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INTRODUCTION 

This  report  describes  the  implementation  of  a  full-wave  fields  program  developed  for 
calculations  at  ELF  (Pappert  and  Shockey,  1977).  The  program  to  be  described  here 
incorporates  modifications  to  the  original  code  for  use  at  VLF,  including  allowance  for  multiple 
modes.  Other  changes  relate  to  improving  compatibility  of  the  basic  program  setup  with  that  of 
other  programs  in  the  Defense  Nuclear  Agency  repertoire. 

The  program  is  written  in  VAX  FORTRAN  for  use  with  Digital  Equipment  Corporation  s 
VAX/VMS  operating  systems.  VAX  FORTRAN  is  based  on  American  National  Standard 
FORTRAN-77  (ANSI  X3.9-1978)  and  includes  support  for  programs  that  conform  to  the 
previous  standard  (ANSI  X3 .9-1 966)  as  well  as  numerous  extensions  to  the  ANSI  standard.  The 
WAVFLD  program  uses  some  extensions  of  the  ANSI  standard,  including  namelist-directed 
input/output  and  additional  data  typing,  such  as  COMPLEX*  16. 
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DESCRIPTION  OF  INPUT 


PROGRAM  CONTROL 

Program  execution  is  controlled  by  a  series  of  character  strings  which  indicate  the  type  of 
data  to  follow.  The  control  strings  must  begin  in  column  one  of  the  record  and  can  be  in 
uppercase  or  lowercase  letters.  The  ionospheric  and  collision  frequency  profiles  can  be  input  by 
either  of  two  different  methods.  If  either  profile  is  strictly  exponential,  it  can  be  read  in  by 
means  of  variables  in  the  NAMELIST  input  signaled  by  the  control  string  NAME.  If  either 
profile  is  non-exponential,  it  must  be  read  in  by  means  of  a  tabular-type  formatted  input  signaled 
by  the  PROFILE  or  COLFREQ  control  strings.  All  of  these  control  strings  are  described  below, 
and  are  summarized  in  Table  1. 


Table  1.  Summary  of  control  strings. 


NAME 

initiates  reading  of  NAMELIST  input 

PROFILE 

initiates  reading  of  ionospheric  profile 

COLFREQ 

initiates  reading  of  ionospheric  collision-frequency  profile 

DATA 

initiates  reading  of  propagation  path  data 

SW  XMT 

provides  same  function  as  DATA 

QUIT 

indicates  end  of  the  input 

NAME  Signals  that  NAMELIST  data  follow.  In  this  program  the  NAMELIST  name 
is  DATUM.  The  NAMELIST  variables  are  described  in  the  next  section. 

PROFILE  Initiates  reading  of  the  ionospheric  charged-particle  profile  data  used  to 
model  the  upper  boundary  of  the  earth-ionosphere  waveguide.  This  allows 
for  using  a  non-exponential  ionospheric  profile.  The  PROFILE  control 
string  is  followed  by  an  alphanumeric  record,  which  is  used  to  identify  the 
profile.  It  can  contain  up  to  80  characters  of  information.  The  profile  is  input 
starting  at  the  top  of  the  ionosphere,  and  the  records  must  be  input  in 
descending  order  of  height.  The  profile  records  contain  the  height,  in 
kilometers,  and  the  species  densities,  in  particles  per  cubic  centimeter,  at  that 
height.  A  format  of  [F7.2,4X,5(1X,E9.2)]  is  used  to  read  the  records.  In  the 
integration  of  the  reflection  elements  through  the  ionosphere,  the  program 
interpolates  exponentially  between  input  values.  The  profile  should  contain 
sufficient  data  to  define  the  ionospheric  structure  with  height.  For  example, 
an  exponential  profile  should  consist  of  only  the  top  and  bottom  heights  and 
densities.  Many  regularly  spaced  heights  tend  to  slow  the  integration.  A 
purely  exponential  ionospheric  profile  (electrons  only)  may  also  be  specified 
in  the  NAMELIST  input  by  the  variables  BETA,  HPRIME,  and  SCLHTS.  A 
maximum  of  101  heights  can  be  specified  and  a  maximum  of  five  species  per 
height  may  be  input.  The  number  of  species  to  be  used  is  determined  from 


the  control  string.  If  the  ninth  column  is  blank,  one  species  is  assumed.  If  a 
value  is  in  this  column,  it  is  used  for  the  number  of  species.  In  the  special  case 
of  three  species,  only  two  species  are  specified  at  each  altitude.  The  first  is 
assumed  to  be  electrons  and  the  second  positive  ions.  The  third  species, 
negative  ions,  is  calculated  in  the  program  by  subtracting  the  electron  density 
from  the  positive-ion  density  to  preserve  charge  neutrality.  All  three  species 
are  listed  in  the  output.  If  the  value  of  any  species  is  less  than  1.0d-20,  it  is  set 
to  that  value.  The  end  of  the  profile  is  indicated  by  a  height  less  than  zero. 

COLFREQ  Initiates  reading  of  the  collision-frequency  profile  via  tabular  input.  This 
allows  for  use  of  a  non-exponential  collision  frequency.  If  the  tabular  input  is 
used,  it  overrides  the  NAMELIST  variables  COEFNU  and  EXPNU. 
Collision  frequencies  for  all  species  must  be  input.  The  ionospheric  profile 
must  precede  the  collision-frequency  profile  because  the  number  of  species  to 
be  used  is  determined  from  the  PROFILE  control  string.  As  with  the 
ionospheric  species  profiles,  the  profile  is  input  starting  at  the  top  of  the 
ionosphere,  the  records  must  be  input  in  descending  order  of  height,  and  the 
program  interpolates  exponentially  between  input  values.  There  need  be  no 
correspondence  between  the  altitudes  used  to  define  the  charged-particle 
profile  and  the  collision-frequency  profile.  Also,  any  species  value  less  than 
I.0d-20  is  replaced  by  that  value,  the  same  format  is  used  to  read  the  records, 
and  the  end  of  the  profile  is  indicated  by  a  height  less  than  zero.  The  profile 
records  contain  the  height  in  kilometers  and  the  collision  frequencies  in 
collisions  per  second.  A  maximum  of  25  heights  can  be  specified,  and  a 
maximum  of  5  species  per  height  may  be  input.  A  strictly  exponential 
collision  frequency  mav  be  specified  in  the  NAMELIST  input  by  the  variables 
COEFNU  and  EXPNU. 

DATA  Signals  that  the  propagation-path  data  follow.  The  format  of  these  data  is  that 

which  is  produced  by  MODEFNDR  (Shellman,  1986)  and  the 
SEGMENTED  WAVEGUIDE  PROGRAM  (Ferguson  and  Snyder,  1987). 
The  first  record  contains  the  data  set  identification.  This  is  followed  by  sets  of 
mode  constants,  one  for  each  path  segment.  The  first  record  for  each  segment 
contains  the  starting  distance,  frequency,  magnetic  azimuth,  codip,  and 
intensity,  and  the  ground  conductivity  and  permittivity.  This  record  is 
followed  by  the  mode-constant  records,  one  pair  for  each  mode,  containing 
the  following  information: 

1  THETA  I  T1  T2 

2  THETA  I  T3  T4, 

where  1  and  2  arc  sequencing  indices  and  THETA  is  the  complex  cigcnangle 
at  the  ground.  For  a  discussion  of  I  and  the  four  T  quantities,  see  Ferguson 
and  Snyder,  1980.  The  list  of  modes  for  each  segment  is  terminated  by  a 


blank  record.  After  reading  the  data  for  each  segment,  the  program  begins 
the  WAVFLD  calculations.  The  list  of  segments  is  terminated  by  a  record 
with  the  starting  distance  set  equal  to  40.0. 

SW  XMT  Provides  the  same  function  as  DATA  control  string. 

QUIT  Indicates  the  end  of  the  input.  The  execution  of  the  program  is  terminated. 

NAMELIST  VARIABLES 

The  NAMELIST  variables  are  described  in  tabular  form  below.  The  initial  or  default 
values  are  given  in  Table  2. 

Table  2.  Namelist  variables  and  initial  values. 


Name 

Value 

Units 

nprof 

0 

beta 

0.0 

1/km 

h  prime 

0.0 

km 

sclhts 

5.0 

h 

0.0 

km 

alpha 

3.1 4c -4 

1/km 

prec 

3.0e-5 

debug 

0 

- 

nelect 

0 

- 

coefnu 

1.816ell,4*4.54c9 

coll/sec 

expnu 

5* -0.15 

1/km 

charge 

-1.0, 1.0, -1.0, 1.0,- 1.0 

mratio 

1.0,4*5.8c4 

topht 

100.0 

km 

Iwstht 

0.0 

km 

itr 

0 

- 

maxitr 

10 

- 

dtheta 

(5.0d-2,1.0d-2) 

degrees 

lub 

(5.0d-2,5.0d-3) 

degrees 

nprint 

1 

- 

savplt 

.true. 

- 

.'1 
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NPROF 


BETA 


HPRIME 


SCLHTS 


ALPHA 


PREC 


DEBUG 


COEFNU 


EXPNU 


A  flag  controlling  usage  of  exponential  electron-density  profile.  Setting  the 
flag  non-zero  indicates  usage  of  exponential  electron-density  profile.  Setting 
the  flag  to  zero  indicates  non-usage  of  exponential  electron-density  profile. 

The  height  variation  of  the  electron  density  for  exponential  ionospheric 
profiles  in  inverse  kilometers. 

The  reference  height  of  the  electron  density  for  exponential  ionospheric 
profiles  in  kilometers. 

The  number  of  scale  heights  above  HPRIME  at  which  the  top  of  the 
exponential  electron-density  profile  is  defined. 

Altitude  in  kilometers  at  which  the  modified  refractive  index  is  unity. 

Earth  curvature  coefficient  in  inverse  kilometers.  ALPHA  is  defined  as 
2/radius  of  the  earth.  For  a  flat  earth,  use  ALPHA  =  0. 

Accuracy  to  be  maintained  locally  in  the  numerical  integrations.  It  is  usually 
taken  to  be  the  default  value  of  3.0e-5. 

A  flag  controlling  the  amount  of  diagnostic  output  desired. 

DEBUG  =  0  suppress  all  printout 
DEBUG  =  1  print  height  gains 

DEBUG  =  2  print  information  on  initial  boundary  conditions  and 
number  of  integration  steps 
DEBUG  =  3  print  extensive  debug  output 

The  collision  frequency,  in  collisions  per  second,  at  the  ground.  It  is  used  with 
EXPNU  to  specify  an  exponential  collision  frequency. 

The  height  variation  in  inverse  kilometers.  It  is  used  to  specify  an  exponential 
collision  frequency. 


CHARGE  The  charge  of  each  species.  For  an  electron  the  charge  is  -1. 

MRATIO  The  ratio  of  the  mass  of  the  species  to  the  mass  of  an  electron.  For  an 
electron,  the  ratio  is  1. 

TOPHT  Greatest  height  in  the  ionosphere  where  field  strengths  are  desired. 

LWSTHT  Lowest  height  in  the  ionosphere  where  field  strengths  arc  desired. 


ITR 


A  flag  to  control  iterations.  If  ITR=  1,  iterations  are  performed  to  refine  the 
input  eigenangle  (THETA). 


MAXITR 


DTHETA 


LUB 


NPRINT 


SAVPLT 


The  maximum  number  of  iterations  allowed  in  finding  a  modal  solution.  The 
iterations  are  terminated  for  the  current  angle  when  MAXITR  is  exceeded. 


The  complex  incremental  change  in  the  modal  solution  for  computing  the 
derivative.  It  is  input  as  a  pair  of  real  numbers. 


A  complex  number  used  to  terminate  the  iterative  process.  The  real  and 
imaginary  parts  of  LUB  are  used  to  test  the  iterative  change  in  the  real  and 
imaginary  parts  of  the  modal  solution.  The  iteration  is  stopped  when  the 
change  is  less  than  or  equal  to  LUB  in  both  real  and  imaginary  parts.  It  is 
input  as  a  pair  of  real  numbers. 


A  flag  controlling  the  amount  of  output. 

NPRINT  =  0  no  output  is  generated. 

NPRINT  =  1  the  total  fields  are  printed. 
NPRINT  =  2  the  individual  modes  are  printed. 


A  logical  flag  for  writing  data  to  logical  unit  10  for  further  processing  and/or 
plotting.  This  variable  allows  the  user  to  submit  batch  or  background  jobs  to 
generate  the  data  to  be  plotted. 

SAVPLT  =  .TRUE.  data  are  written 
SAVPLT  =  .FALSE.  no  data  are  written 


>1 


>1 


% 


SAMPLE  PROBLEM 


A  sample  case  will  be  discussed.  The  data  for  this  case  were  generated  using 
MODEFNDR  (Shellman,  1986).  The  printed  output  for  this  MODEFNDR  run  is  shown  in 
Figure  1.  This  figure  is  included  here  so  the  user  may  see  the  correspondence  between  inputs  to 
MODEFNDR  and  the  inputs  to  WAVFLD  (shown  in  Figure  2).  The  profile  specifications  are 
identical  in  the  two  programs.  The  data  lines  between  ‘DATA’  and  the  blank  line  just  before  ‘R 
40’  are  all  output  from  MODEFNDR.  Immediately  after  the  DATA  control  string  is  the  data 
identification  string  which  followed  the  ID  string  in  Figure  1.  The  frequency,  geomagnetic  field 
parameters,  and  ground  conductivity  parameters  are  encoded  in  the  line  after  the  identification 
string.  This  is  followed  by  pairs  of  lines  containing  the  mode  parameters. 

The  printed  output  from  WAVFLD  is  shown  in  Figure  3.  This  output  shows  the 
ionospheric  profile  data  and  the  values  of  the  NAMELIST  variables.  After  the  DATA  control 
string,  the  parameters  of  the  individual  modes  are  shown.  Initially,  the  basic  parameters  of 
attenuation  rate,  normalized  phase  velocity,  and  excitation  factor  are  printed.  The  number  of 
integration  steps  for  one  integration  and  the  resulting  value  of  the  complex  mode  equation  are 
printed.  The  sample  case  calls  for  iteration  of  the  input  eigenangles,  so  the  integrations  continue 
until  the  incremental  change  in  the  eigenangle  is  less  than  LUB.  For  the  first  mode,  this  happens 
after  two  iterations.  In  the  same  fashion  the  remaining  modes  are  processed.  The  electric  and 
magnetic  components  for  each  mode  are  adjusted  to  correspond  to  the  calculated  excitation 
factor  and  the  total  fields  are  computed  by  summation.  The  resulting  values  of  the  fields  as  a 
function  of  altitude  are  printed.  In  addition,  the  relative  Joule  heating  in  watts  per  unit  volume 
and  the  components  of  the  Poynting  vector  are  printed  under  the  headings  ‘q’,  ‘sx\  ‘sy’  and  ‘sz’. 
The  angles  that  the  Poynting  vector  makes  with  respect  to  the  vertical  are  shown  under  the 
headings  ‘alpha’  and  ‘beta.’  Representative  plots  for  these  test  data  are  shown  in  Figures  4 
through  8.  It  should  be  noted  that  the  user  may  obtain  other  parameters  (such  as  described  by 
Pappert  and  Shockey,  1977)  by  modifications  to  the  main  routine. 
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PROFILE 
TEST  CASE 


120.00 

5.80E 

02 

112.00 

1.10E 

03 

110.00 

1.30E 

03 

106.00 

1.70E 

03 

104.00 

1.90E 

03 

102.00 

1.98E 

03 

100.00 

2.00E 

03 

99.00 

1.95E 

03 

0.00 

1.83E- 

■12 

-99. 

COLFREQ 

120.00 

1.00E 

04 

104.00 

3.00E 

04 

0.00 

1.82E 

11 

WAVFLD  test  data 
NAME 
&DATUM 
FREQ= 1 7.8, 

AZIM=87 .801 ,C0DIP=6. 593 ,MAGFLO= .52220-4, 
SIGMA-1. 0D-5.EPSR-5.0, 
RANGER=80,90,RANGEI=0, -5, 

SEND 

THE  TOP  OF  THE  PROFILE  IS  SET  TO  BE-- 
97.11  1.01E+03 

AT  THE  TOP  OF  THE  PROFILE  B  =  3.519E+01 

THE  BOTTOM  PROFILE  HEIGHT  IS  =  67.047 

OMEGA-R  EQUALS  2.5E5  AT  HT  =  87.06 


OMEGA-R  = 
83.696 
74.409 
74.409 


2.503E+05  AT  HOFWR  = 


-0.665 

-7.662 

-7.662 


84.804 


3 . 519E+01  OMEGA-R  =  3.806E+07 

67.047  AT  B(CUTOFF)  =  1.000E-04 

87.06 
87.06 
-0.101 


74.409  -7.662  USED 


74.409  -7.662  DELETED 


R  0.00  F  17, 

.8000  A 

87.801 

C  6.593 

M  0.522E- 

04  S  1 . 000E-05  E  5. 

,0  T  87. 

MODE  THETA 

ATTEN 

VOVERC 

WAIT  MAG 

WAIT  ANG  THETAP 

1  84.804 

-0.101 

0.523 

0.99624 

-64.578 

2.686  89.892 

-4.915 

2  83.696 

-0.665 

4.164 

0.99812 

-23.736 

2.850  88.836 

-3.622 

3  74.409 

-7.662 

117.736 

1.02091 

16.911 

-1.104  75.698 

-8.396 

QUIT 

***THE  CALCULATIONS  ARE  COMPLETE*** 

FORTRAN  STOP 


Figure  1.  Sample  output  from  MODEFNDR.  The  data  obtained  in  this  case 
are  used  in  the  sample  run  for  WAVFLD. 
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PROFILE 

TEST  CASE 

120.00 

5.80E 

02 

112.00 

1.10E 

03 

110.00 

1.30E 

03 

106.00 

1.70E 

03 

104.00 

1.90E 

03 

102.00 

1.98E 

03 

100.00 

2.00E 

03 

99.00 

1.95E 

03 

0.00 

1.83E- 

■12 

-99. 

COLFREQ 

120.00 

1.00E 

04 

104.00 

3.00E 

04 

0.00 

1.82E 

11 

-99. 

NAME 

&DATUM 

T0PHT=120.0,LWSTHT=0.0, 

ITR-1, 

SEND 

DATA 

WAVFLD  test  data 

R  0.000  F  17.8000  A  87.801  C  6.593  M  0.522E-04  S  1.000E-05  E  5.0  T  87.1 

1  89.89229  -4.914532  9 . 641 1 4224E - 06  1 . 9681 5890E - 05 - 6 . 95767994 E - 06  2 . 03831610E-06 

2  89.89229  -4.914532  5.02810508E-06-1 . 14897985E-05  1 . 0101 0334E+00  5 . 58272935E-03 

1  88.83618  -3.622151  6 . 90992048E-04  2 . 3 2309 1 73 E - 03 - 2 . 42370220E - 07 - 1 . 95802787E-08 

2  88.83618  -3.622151-1 . 5 1 57826 1 E-05  1 . 88085523 E - 05  1 . 01016688E+00  5 . 75073995E-03 

1  75.69807  -8. 396441-2. 33791590E-01 - 1 .40563220E-01 -5 . 267631 12E- 1 1  4 . 32070248E- 1 1 

2  75.69807  -8. 396441-4. 27416853E-06  3 .35867099E-07  1 . 01005256E+00  8. 58251 192E-03 

R  40. 


Figure  2.  Sample  input  to  the  WAVFLD  program. 
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APPENDIX  A 


SOURCE  LISTING  FOR  PROGRAM  WAVFLD 


wavfld 


name! ist:  nprof,beta,hprime,sclhts, 
h, alpha, prec.de bug, nelect, 
coefnu,expnu,charge,mratio, 
topht,lwstht,itr,maxitr,dtheta,lub, 
npr i nt , savpl t 

nprint=0  gives  no  print 

1  print  of  total  fields  etc. 

2  print  of  individual  modes 

savplt=t  writes  data  to  unit  10  for  processing  by  wfplts 
=f  no  output  to  unit  10 

implicit  real*8  (a-h,o-z) 

common/wfinpt/thetap,freq,azim,codip,magfld,coefnu(5),expnu(5), 

$  topht,lwstht,delht,h, alpha, si gma,epsr,nelect 

common/wf  fl ag/prec , i so , debug 

common/wf  prof/'htl  ist(101)  ,  1  n  1  i  s  t  ( 1 0 1  ,5)  ,  hi  i  st  (25)  ,  cf  1  i  s  t  (  25 , 5  )  , 
$  nrlht,lht,nrmht,mht,charge(5),mratio(5),nrspec 

common/ i trat/dtheta(2),lub(2),maxitr,itr 
common/m  mtx/mll,m21,m31,ml2,m22,m32,ml3,m23,m33 
dimension  en(5),nu(5),label(2),plotx(101,16),ploty(101), 

$  ex(101),  ey(101),  ez(101),  hx(101),  hy(101),  hz(101), 

S  sex(lOl) ,sey(101) ,sez(101) ,shx(101) , shy(lOl)  ,shz(101) , 

S  smll(101),sm21(101),sm31(101), 

S  sml2 ( 101 ) , sm22 ( 101 ) , sm32 ( 101 ) , 

$  sml3(101) , sm23 ( 101 ) , sm33 ( 101 ) 

comp lex* 16  mi/ (0. dO, -l.dO)/, theta, thetap,st,stp,exc, tempi, temp 2, 
$  ex,ey,ez,hx,hy,hz,sex,sey,sez,shx,shy,shz, 

S  smll,sm21,sm31,sml2,sm22,sm32,sml3,sm23,sm33, 

$  mll,m21,m31,ml2,m22,m32,ml3,m23,m33,el,e2,e3,hl,h2,h3 

real*8  magfld,nu,mratio,lnlist,lwstht,lub 
real*4  plotx,ploty,capy,frq,rho,azm,cdp.sig,tpht,lwht,attn 
logical  f i rst/ . true ./, savpl t 
integer  debug, hdate(3) , h t i me ( 2 ) 


character*8  cti me, branch, bed ( 1 0 ) , label/’ 
character*9  edate 


, ’ total 


equivalence  (cdate,hdate),(ctime,htime) 
namel ist/datum/ h, alpha, prec, debug, nelect, 

S  coefnu,expnu,charge,mratio,nprof,beta,hprime,sclhts, 

$  topht,lwstht,itr,maxitr,dtheta,lub, 

$  npr i nt , savpl t 

data  dtr/1 .745329252d-2/,npflag,nuflag/2*0/, 

S  n prof/0/ , be ta , hpr i me/ 2*0 . dO/ , sc  1 ht s/5 . dO/ , 

S  nprint/1/, savpl t/ . true./ 

data  freq/0.d0/,azim/0.d0/,codip/0.d0/,magfld/0.d0/, 

S  coefnu/ 1 . 81 6d 1 1 , 4*4 . 54d9/ , expnu/5*- . 1 5d0/ , 

S  t oph t/ 1 00./, lwstht/0. dO/, nrht/101/,h/0.d0/, alpha/3. 14d - 4/ 

S  sigma/5. d0/,epsr/ 81. dO/, nelect/ 0/, 

$  prec/3 . d- 5/ , debug/0/ , 

$  Charge/-l.d0,l.d0,-l.d0,l.d0,-l.d0/,mratio/l.d0,4*5.8d4/, 

$  itr/0/,maxitr/10/,dtheta/5.d-2,l.d-2/,lub/5.d-2,5.d-3/ 


*,rf,**i 


call  date(cdate) 
call  time(ctime) 
print  1000,cdate,ctime 


read(5, 180, end 
print  181, bed 
if(bcd(l)  .eq. 

bcd(l)  .eq. 
i f ( bed ( 1 )  .eq. 

bcd(l)  .eq. 
i f ( bed ( 1 )  .eq. 

bcd(l)  .eq. 
if(bcd(l)  .eq. 

bcd(l)  .eq. 
i f ( bed ( 1 )  .eq. 

bcd(l)  .eq. 
i f (bed ( 1 )  . eq . 

bcd(l)  .eq. 
i f ( bed ( 1 )  .eq. 

bcd(l)  .eq. 
go  to  910 


=999)  bed 

’quit  ' 

’QUIT  ’) 

’name 

’NAME  ’) 

’profile  ’ 
’PROFILE  ’) 
’colfreq  ’ 

’ COLFREQ  ’) 
’data 

’DATA  ’) 

’ipsq  xmt’ 

’  IPSQ  XMT’ ) 
’sw  xmt’ 
’SW  XMT’) 


go  to  999 
.  or . 

go  to  20 
.  or . 

go  to  30 
.  or . 

go  to  50 
.or. 

go  to  60 
.  or . 

go  to  60 
.  or. 

go  to  60 


namelist  input 

read(5, datum) 

print  datum 

capk= 1 . - . 5*al pha*h 

i f (nprof  .eq.  0)  go  to  25 

if (beta*hprime  .eq.  O.dO)  go  to  913 

npflag  =  l 

nn=hprime+sclhts/beta+. 5d0 
nrspec=l 

c0=dlog(7.85535d-05*coefnu(l))-beta*hprime 
cl=beta+expnu ( 1 ) 
nrl ht=2 

htl  i st ( 1 ) =topht 
htl  i st ( 2 ) =0 . dO 
do  24  1=1,2 

1  n  1  i  st ( 1 , 1 )  =c0  +  c l*ht 1  i  st  ( 1  ) 

number  of  points  set  to  101 
delht=(topht-lwstht)/100.d0 

if(nuflag  .eq.  1)  go  to  10 
nrmht=2 

hi  i  st ( 1 ) = topht 

hi  i st ( 2 ) =0 . dO 

do  27  k=l,nrspec 

cflist(2,k)=dlog(coefnu(k)) 

cflist(l,k)=cflist(2,k)+topht*expnu(k) 

go  to  10 

profile  input 

decode(9, 170, bed)  branch, nn 
npfl ag  =  1 
nr spec=nn 

if(nrspec  .eq.  0)  nrspec=l 
read(5,180)  bed 
print  181, bed 
1=0 


31  read(5, 140)  ht,en 

i f (ht  .It.  0  .d0)  go  to  39 
i f (1  . eq .  101 )  go  to  912 
i  f  ( 1  .eq.  0)  go  to  32 
if(ht  .ge.  htlist(l))  go  to  911 

32  1=1+1 

if(nrspec  .eq.  3)  en (3 ) =en ( 2 ) - en ( 1 ) 
print  141 , ht ,( en (k) , k=l , nrspec ) 
htl i st (1 )=ht 
do  33  k= 1 , nr  spec 

33  1 n 1 i s  t ( 1 ,k)=dlog(dmaxl(en(k),l.d-20)) 
go  to  31 

39  nrlht=l 

go  to  10 

c 

c  collision  frequency  profiles 

50  nuflag  =  1 
1=0 

51  read(5,140)  ht,nu 

i  f  (ht  .It.  0  .d0)  go  to  59 
i f ( 1  . eq  .  25 )  go  to  912 
i  f  ( 1  .eq.  0)  go  to  52 
if(ht  .ge.  hlist(l))  go  to  911 

52  1=1+1 

print  141,  ht ,  ( n u ( k ) , k= 1 , nrspec) 
hi i st ( 1  )  =ht 
do  53  k  =  1 , nrspec 

53  c  f 1 i st ( 1 ,k)=dlog(dmaxl(nu(k),l.d-20)) 
go  to  51 

59  nrmht=l 
go  to  10 

c 

c  wavefields  calculations 

60  i f ( n p f 1  a g  .eq.  0)  go  to  915 

read(5,  180)  bed 

61  read(5,1010)  r  ,  f , a , c , b , s , e 
if(r  .eq.  40. d 0 )  go  to  300 
i  f  (  s  .  eq  .  0  .  dO )  go  to  61 
rho  =  r 

freq=f 
az i m=a 
cod i p=c 
magf 1 d=b 
b=b* 1 . d4 

d  c 1 =  dsin(codip*dtr)*dcos(azim*dtr) 
dcm=  dsin(codip*dtr)*dsi n(azim*dtr ) 
dcn=-dcos(codip*dtr) 
s  i  gma  =  s 
epsr  =  e 

omega=6. 283 1 85307 18d03*freq 

coefx=3.182357d09/omega**2 

capy= 1 . 7  58  7  96d 1 1 *magf 1 d/ omega 

wn=20 . 95845d0*freq 

const=- freq 

aconst=-8. 686d0*wn 

econst=- . 035d0*wn 

nsum=  1 

nm=0 


A-4 


V.] 


if(nprint  .le.  1)  go  to  62 
print  160 
print  181, bed 

print  101 1 , r , f , a , c , b , s  , e , h 

read (5, 1023)  indxi,trl,til,tmprl,tmpil,tmpr2,tmpi2 
if(trl  .eq.  O.dO)  go  to  91 
nm=nm+l 

re  ad  (5, 1023)  indx2,tr2, ti2,tmpr3,tmpi3,tmpr4,tmpi4 
if(trl  .ne.  tr2  .or.  til  .ne.  t i 2 )  go  to  916 
thetap=dcmpl x ( tr 1 , t i 1 ) 
if(tmprl  .eq.  O.dO  .and.  tmpil  .eq.  O.dO)  tmprl=l.d-20 
templ=dcmpl x (tmprl , tmpi  1 ) 
i  f ( i n  d  x 1  .eq.  1)  go  to  63 

if(tmpr3  .eq.  O.dO  .and.  tmpi3  .eq.  O.dO)  tmpr3=l.< 
templ=dcmpl x ( tmpr3 , tmp i 3 ) 


O.dO)  tmpr3= 1 . d - 20 


stp=cdsin(thetap*dtr) 
sr =stp 
s i =stp*mi 
atten=aconst*si 
voverc=l . dO/sr 

exc=(-2.124292957d0,0.d0)*templ*stp*stp 

excr=dcmplx(0.d0,econst)*exc 

exci=dcmplx(0.dO,econst)*exc*mi 

wm=10.d0*dlogl0(excr**2+exci**2) 

wa=datan2 (exci ,excr) 

if(h  .eq.  O.dO)  go  to  64 

st=stp*capk 

theta=mi *cdl  og(cdsqrt(l.dO-st*st)-mi*st)/dtr 

go  to  65 

theta=thetap 

i f (npr  i  nt  . 1 e .  1 )  go  to  66 
print  160 
print  181, bed 

print  1011 ,r,f ,a,c,b,s,e,h 

print  1040, nm, theta, at ten, vo verc,wm,wa,thetap 

call  wavfl d(ex, ey, ez , hx, by, hz) 

lht-0 

ht=topht 

temp2=exc 

i f ( npr i nt  .gt.  1)  go  to  75 

if(nprint  .gt.  0  .and.  nsum  .eq.  1)  go  to  80 

print  900, label (nsum) 

do  90  k=l,nrht 

j=nrht+l -k 

if(nsum  .eq.  1)  go  to  81 

el=sex(j) 

e2=sey(j) 

e3=sez( j) 

h  1  =  s  h  x  ( j  ) 

h2  =  shy ( j ) 

h3  =  shz ( j ) 

go  to  82 

el«ex(j) 

e2=ey( j ) 

e  3  =  e  z  ( j ) 


h 1  =  h  x ( j ) 
h2=hy(j) 
h3  =  hz (j ) 

if(nm  .eq.  1)  go  to  83 
mll=smll ( j) 
m21=sm21 ( j ) 
m31=sm31 ( j ) 
ml2=sml2( j ) 
m22=sm22 ( j ) 
m32=sm32 ( j ) 
ml3=sml3 ( j ) 
m,23  =  sm23  ( j ) 
m33  =  sm33  ( j ) 

if(nsum  .eq.  2)  go  to  84 

sex( j ) =sex( j )+ex(j )*temp2 

sey( j)=sey(j)+ey(j)*temp2 

sez( j)=sez(j )+ez( j)*temp2 

shx(j)=shx(j)+hx(j)*temp2 

shy( j)=shy(j)+hy(j )*temp2 

shz( j)=shz(j)+hz(j)*temp2 

go  to  84 

call  tmtrx(ht) 

sex( j)=ex( j)*temp2 

sey( j)=ey( j)*temp2 

sez( j)=ez( j)*temp2 

shx(  j)=hx(.j)*temp2 

shy( j)=hy( j)*temp2 

shz(j)=hz(j)*temp2 

sml 1 ( j ) =  m 1 1 

sm21 ( j )=m21 

sm31 ( j ) =m3 1 

sml2(j)=ml2 

sm22 ( j )=m22 

sm32 ( j ) =m32 

sml3(j)=ml3 

sm23 ( j  j  =m23 

sm33 ( j ) =m33 

q  and  s  are  relative 

q=dcmpl x(0.d0, const) *(dconjg (ml I*el+ml2*e2+ml3*e3)*el 
$  +dcon jg (m21*e l+m22*e2+m23*e3 ) *e2 

$  +dconjg(m31*el+m32*e2+m33*e3)*e3) 

suma=suma+q*del ht 
sx=e2*dconjg(h3)-dconjg(h2)*e3 
sy=e3*dconjg(hl)-dconjg(h3)*el 
sz=el*dconjg(h2)-dconjg(hl)*e2 
exm=cdabs (el ) 
eym=cdabs ( e2 ) 
ezm=cdabs (e3 ) 
hxm=cdabs ( h 1 ) 
hym=cdabs ( h2 ) 
hzm=cdabs ( h3 ) 
exa=cdang(el ) 
eya=cdang ( e2 ) 
eza=cdang(e3) 
hxa=cdang(hl ) 
hya=cdang ( h2 ) 
hza=cdang ( h3 ) 
i f (nsum  . eq .  1 )  go  to  87 


ploty{k)=ht 

c  joule  heating  term 

plotx(k,l)=dlogl0(dmaxl(l.d-10,dabs(q  ))) 
c  attenuation 

plotx(k,2)=suma 

c  poynting  vector  magnitudes 

plotx(k,3)=dlogl0(dmaxl(l.d-10,dabs(sx))) 
plotx(k,4)=dlogl0(dmaxl(l.d-10,dabs(sy))) 
pi  otx( k, 5) =dl og 10 (dmaxl ( 1 . d- 10 ,dabs ( sz ) ) ) 
angl=datan2(sz,sx)/dtr 
ang2=datan2(sz,sy)/dtr 

ang3=dacos((dcl *sx+dcm*sy+dcn*sz)/sqrt (sx**2+sy**2+sz**2))/dtr 
c  angles  between  poynting  vectors 

plotx(k,6)=angl 
plotx(k,7)=ang2 
plotx(k,8)=ang3 
call  wfdens(ht,en,nu) 
capx=coefx*en ( 1 ) 
capz=nu ( 1 )/omega 
c  ionospheric  parameters 

pi otx ( k,  9)=dlogl0{dmaxl(l.d-10,dabs(capx))) 
pi otx( k, 1 0 ) =dl og 10 (dmaxl ( 1 .d- 10, dabs (capz) ) ) 
c  magnitude  of  electric  field  vectors 

plotx(k,ll)=dlogl0(dmaxl(l.d-10,exm)) 
plotx(k,12)=dlogl0(dmaxl(l.d-10, eym) ) 
pi otx(k, 13 ) =dl oglO (dmaxl ( 1 .d- 10, ezm) ) 
c  magnitude  of  magnetic  field  vectors 

pi otx( k, 14)=dl og 10 (dmaxl (1 .d-10, hxm) ) 
plotx(k,15)=dlogl0(dmaxl(l . d- 10, hym) ) 
plotx(k,16)=dlogl0(dmaxl(l.d-10,hzm)) 
if(nprint  .gt.  0)  go  to  88 
go  to  89 

87  if(nprint  .le.  1)  go  to  89 

88  print  90 1 , ht , exm, exa , eym, eya , ezm , eza , q , sx , sy , sz , ang 1 , ang2 , 

$  hxm,hxa,hym,hya,hzm,hza 

89  ht=ht-delht 

if(ht  .It.  htl i st (1 ht+1 ) )  lht=lht+l 
i f ( h t  .It.  hi i st (mht  +  1 ) )  mht=mht+l 

90  continue 
c 

if(nsum  .ne.  2)  go  to  62 
go  to  92 

91  nsum=2 
suma*-0 . 

if(nprint  .eq.  0)  go  to  74 
print  160 
print  181, bed 

print  101 1 , r , f , a , c , b , s , e , h 
print  1012,capy,nm 
go  to  74 
c 

92  if(savplt)  then 
c 

if(first)  then 
f i rst= . fal se . 

open(unit=10,type=’new’ , form= ’unformatted’) 
end  if 
rho=r 


A- 7 


f 

W 


frq=freq 

azm=azim 

cdp=codi p 

sig=sigma 

tpht=topht 

1 wh  t  =  1 wstht 

attn=suma 

:  save  plot  data  on  a  file  for  later 

write(lO)  bcd,hdate,htime,frq,rho,azm,cdp,sig,capy, 

$  nrht,tpht,lwht,attn,plotx,ploty 

end  i  f 
go  to  61 

iOO  if(nprint  .gt.  0)  print  160 
go  to  10 

310  print  9910 

go  to  999 

311  print  9911 

go  to  999 

312  print  9912 

go  to  999 

313  print  9913 

go  to  999 

314  print  9914 

go  to  999 

315  print  9915 

go  to  999 

316  print  9916 

399  stop 

140  format(f7.2,4x,5(lx,e9.2)) 

141  format(f8.2,4x,lp5e9.2) 

160  format(lhl) 

161  format(lh  ) 

170  format ( 1 a8 , 1 i 1 ) 

171  format ( lx, la8) 

180  format(10a8) 

181  format ( lx, 10a8) 

900  format (’ 0 ’, a5, ’  Field  strengths  computed  by  WAVFLD’/ 

$  ’0  ht’,6x,’mag  --  x  --  ang  mag  --  y  --  ang’, 

$  3x,’mag  --  z  --  ang  q ’ , lOx , ’ sx ’ , 9x , ’ sy ’ , 9x , ’ sz ’ , 


$  ’0  ht’,6x,’mag  --  x  --  ang  mag  --  y  --  ang’, 

$  3x,’mag  --  z  --  ang  q ’ , lOx , ’ sx ’ , 9x , ’ sy ’ , 9x , ’ sz ’  , 

$  1 2 x , ’alpha’  ,6x, ’beta’ ) 

format ( lx, f 5 . 1 , ’  e ’ , 3 ( 1  pel  1 . 3 , Opf 7 . 3) , 1 p4e 1 1 . 3 , 0p2f 1 1 . 3/ 

$  8x ,  ’h’ , 3 ( 1 p  e 1 1 .3,0pf7.3) ) 

format(’  Additional  plot  identification:  ’,a9,2x,a8/) 
format(lx,f7.3,2x,f8.4,2x,f8.3,2x,f8.3,2x,el0.3,2x,el0.3,2x,f5.1) 
format ( ’ Or ’ , f7 . 3 , ’  f’,f8.4,’  a’,f8.3,’  c’,f8.3,’  m’,f6.3, 

$  ’  s  ’ , IpelO . 3, ’  e ’ , 0 p f 5 . 1 , ’  h ’ , f 5 . 1 ) 

format (’+’, 72x , ’  mag(y) ’ , lpel0.3, ’  modes’, i2) 
format ( i 1 , 2f9. 0, lx,4el5.0) 

format (’ Oi nput  for  nm  =  ’ ,  i 2 ,  ’  :  ’ , i 1 , Op 2 f 10 . 5 , 1 p 2 e 1 6 . 8/ 

S  2  1  x  ,  i 1 , 0 p 2 f 1 0 . 5 , 1 p2el 6 . 8) 

format (’ Omode  real  imag  atten  v/c’,9x, 

S  ’waif’s  exc  real’’  imag’’’/ 

$  lx, i 3 , f 1 1 .5,fl0.5,f9.3,f9.5,fl0.3,f7.3,fl0.5,fl0.5) 

format(/’  Error  in  control  string’) 
format (/’  Heights  in  profile  out  of  order’) 


i£ 


>  *  *  »  ♦»>  »>*  A«.WA*  j*  m*M  tj  ij  i 


9912 

9913 

9914 

9915 

9916 


format (/ ’ 
format (/’ 
format (/ ’ 
format (/ ’ 
format (/  ’ 
end 


Too  many  heights  in  profile’) 

BETA  or  HPRIME  are  0.0’) 

Number  of  output  heights  is  greater  than  101’) 
No  profile  specified’) 

Input  data  out  of  order’) 


r 


$ 

•S' 

•» 


M 
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subroutine  wavfld(ex,ey,ez,hx,hy,hz) 


wavfld  calls  for  the  downward  integration,  and  then  performs  the 
back  substitution  of  normalizing  values  (saved  as  data  by  wfstor). 
field  strengths  are  computed  at  heights  from  topht  to  lwstht  at 
delht  increments  and  are  returned  in  ex,  ey,  ez,  hx,  hy,  hz. 

implicit  real*8  (a-h,o-z) 

common/wf  inpt/theta, f req, azmuth , codip,magfld,ceffnu(5),expnu(5), 

$  topht, lwstht, delht, h, alpha, sigma, epslon,nglect 

common/wf  f 1 ag/precsn , i so , i dbg 

common/wf  save/p (4 , 2) ,m31 ,m32 , m3 3 , ortho , anorm, bnorm, ht , 1 evl 
common/cs/c , s , ci , s  i 
real*8  magfld, lwstht 

compl ex* 16  theta, ex ( 1 ) , ey ( 1 ) , ez ( 1 ) , hx ( 1 ) , hy ( 1 ) , hz ( 1 ) , 

S  p,m31,m32,m33,ortho,c,s,ci,si,b(2),w(4),osum 

jht=topht/delht+1.01d0 
test=(jht-l) *del ht-topht 
if (dabs(test)  .gt.  l.d-4)  go  to  800 
mht=lwstht/delht+1.01d0 
test=(mht-l)*delht-lwstht 
if (dabs(test)  .gt.  l.d-4)  go  to  800 
jht=min0(jht-mht+l , 101) 
mht  =  1 

iteration  to  satisfy  modal  equation 
call  itrate 

combine  solutions  at  ground  so  that  they  satisfy  boundary  condition, 
call  wf  bndy(b) 

perform  back  substitution  of  normalizing  values. 

20  o  sum=0.0 
proda=  1 . 0 
prodb=l .  0 
i  h  t = m  h  t 
call  wf  load 
go  to  25 

21  o  sum=o  sum*anorm/bnorm+ortho 
proda=proda*anorm 

if(proda  .It.  1.0d-30)  proda=0.0 
p rod b= prod b* bnorm 
call  wf  load 
do  23  j-1,4 

p(j,2)=(p(j,2) -o  sum*p ( j , 1 ) )*prodb 
23  p(j, l)=p(j,l)*proda 

compute  field  strengths  at  profile  heights. 

25  do  26  j-1,4 

26  w(j)=p(j,l)*b(l)+p(j,2)*b(2) 
ex ( i  ht )=w( 1 ) 

ey ( i ht)=-w(2) 

ez ( i ht )= - ( s*w( 4) +m3 l*w( 1 ) -m3 2*w( 2 ) )/ ( 1 ,0+m33) 
hx( iht)=w(3) 
hy ( i ht ) =w( 4) 
hz ( i ht ) =- s*w(2) 


I 
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I 
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iht=iht+l 

i f ( i h t  . le.  jht)  go  to  21 
if(levl  .ne.  0)  print  903, lev! 
return 

800  print  902 
stop 

902  format (’ Oerror  in  wavfld’/ 

$  ’  delht  does  not  divide  topht-lwstht  evenly’) 

903  format (’ 01 evl  not  zero:  levl=’,i3) 


subroutine  wf  i ntg (topht , 1 wstht , del ht , if 1 ag) 

wf  intg  performs  the  integration  of  the  p  matrix  down  through  the 
ionosphere,  using  the  techniques  given  by  pitteway. 
accuracy  is  maintained  by  adjusting  the  stepsize  so  that  the 
p  matrix  is  computed  with  sufficient  accuracy. 

i f 1 ag =0  integ  for  theta  only 
i f 1 ag  =  1  integ  for  theta  and  theta-dtheta 

implicit  real*8  (a-h,o-z) 
common/wf  fl ag/precsn , i so, idbg 
common/p  mtx/p(16) ,pi (16) 

commcn/wf  save/p  save(16),m31  sav,m32  sav,m33  sav, 

$  ortho, anorm, bnorm, ht ,  1  evl 

common/m  mtx/m(3,3) 

common/wf  prof/enht(101),enlog(101,5),collht(25),collfr(25,5), 

$  nrl ht , 1 ht , nrmht , mht,charge(5),ratiom(5),nrspec 

integer  svflag 
real*8  lwstht 

complex*16  m31  sav,m32  sav,m33  sav, ortho, m 

dimension  prevp ( 1 6 ) , tempp ( 16 ) , dpdh ( 1 6 ) , p v  dpdh ( 16) ,dpidh(16) 

minimum  step-size  allowed 

data  epsht/5 .d -4/ , dhmi n/1 . d - 3/ 

call  i ni t  t 

call  t  mtrx(top  ht) 

call  wf  init(p) 

call  p  deriv(p,dpdh) 

if(iflag  .eq.  0)  go  to  11 

call  ti  mtrx 

call  wf  init(pi ) 

cal  1  p  deriv(pi ,dpidh) 

11  continue 

i steps=0 
kmax=0 
1 evl =0 
ht=topht 

cal  1  xfer(p,p  save, 16) 
m31  sav=m(3,l) 
m3 2  sav=m(3,2) 
m33  sav=m(3,3) 
call  wf  stor 
wf ht=topht -del ht 
del h  2  =  0 . 1 2  5d0*del h t 
svf 1 ag=0 

determine  next  stepsize  to  use. 
if(svflag  .eq.  1)  delh2=savdh2 
svf 1 ag=0 
nodbl =0 
h  1 0  =  h  t 

cal  1  xfer(p, prevp, 16) 
call  xfer(dpdh,pv  dpdh, 16) 
htl im=wfht 

if (enh t ( 1 ht+ 1 )  .gt.  htlim+epsht)  htl im=enht ( 1  ht+1 ) 


if (col lht(mht+l)  .gt.  htl im+epsht)  htl im=col 1 ht (mht+1 ) 

if (htO-del h2 .ge.htl im+epsht)  go  to  50 

savdh2=del h2 

svfl  ag  =  l 

del h2=ht0-htl im 

perform  next  integration  step. 

50  call  wf  step(p,dpdh,ht,delh2,0) 
cal  1  xfer(p, tempp, 16) 
m3 1  sav=m( 3 , 1 ) 
m32  sav=m(3,2) 
m33  sav=m(3,3) 
h  t  =  h  1 0 

call  xfer(prevp, p, 16) 
call  xfer(pv  dpdh , dpdh , 1 6 ) 
del h=0 . 5  *  d  e 1 h2 

call  wf  step ( p , dpdh , ht , del h , 1 ) 
cal  1  p  deri v(p, dpdh) 
call  wf  step (p , dpdh , ht , del h , 2) 
check  accuracy  of  result. 
pmax=0 . 0 
do  85  j  =  l ,  16 
pabs=dabs(p( j ) - tempp ( j ) ) 
if(pmax  .It.  pabs)  pmax=pabs 
85  continue 

adjust  stepsize  if  necessary. 
if(pmax  .It.  precsn)  go  to  100 

i f (del h  .gt.  dhmin)  go  to  95 
if(kmax  .eq.  0)  print  900,  ht 
kmax=l 
go  to  100 
95  continue 

del h  2  =  0 . 5*delh2 
nodbl = 1 

if(pmax  .It.  1 0 . 0*prec sn )  go  to  99 
del h2=0 . 25*del h2 
nodbl =0 
99  continue 
ht-htO 

cal  1  xfer ( prevp , p , 16) 
call  xfer(pv  dpdh , dpdh , 1 6) 
svfl ag=0 
go  to  50 

100  cal  1  wf  seal (p,0) 

call  xfer(p,p  save, 16) 

if(ht  .It.  wfht+epsht)  call  wf  stor 

cal  1  p  deri v(p, dpdh) 

if(iflag  .eq.  0)  go  to  72 

h  t  =  h  1 0 

call  wf  step(pi ,dpidh,ht,delh,3) 

call  p  deriv(pi , d p i d h ) 

call  wf  step(pi , dp idh , ht , del h , 4 ) 

cal  1  wf  seal (pi , 1 ) 

call  p  deriv(pi , d p i d h ) 

72  continue 


i  s  t  e  p  s  =  i steps  +  1 
i f ( i d bg  .eq.  0)  go  to  73 
idi v=i steps/50 

if(isteps  .eq.  50*idiv)  print  902 , i steps , ht 
73  continue 

if(no  d  b 1  .eq.  0  .and.  pmax  .It.  0.1*precsn)  del h  2  =  2 . 0*del h2 

check  integration  and  profile  heights. 
if(ht  .It.  1 wstht+epsht)  go  to  80 
if(ht  .It.  wfht  +  epsht)  wf ht=wfht - del ht 
if(ht  .It.  enht ( 1 h t+ 1 ) +epsht )  1 ht=l ht+1 

if(ht  .It.  col  1 ht (mht+1 ) +epsht )  mht=mht+l 
go  to  10 

80  print  901,isteps 
return 

900  format (’ Omi n imum  stepsize  used  at  h t  = ’ , 1 pe 1 4 . 5) 

901  format (’ 0 i 4 , ’  integration  steps  used  in  wavfld’) 

902  format( ’0’ , i 4 , ’  integration  steps,  ht=’,f9.4) 
end 


subroutine  wf  seal (pp, if lag) 

wfscal  scales  and  orthogonal i zes  the  solution  vectors  p. 

this  scaling  must  later  be  removed  to  yield  correct  (unsealed) 

solutions. 

implicit  real*8  (a-h,o-z) 

common/wf  save/p  save(16),m31  sav,m32  sav,m33  sav, 

$  o  sum, aprod, bprod, ht , 1 evl 

common/save/p  etc(27,101) 

complex*16  p(4,2),m31  sav,m32  sav,m33  sav,o  sum, ortho 
dimension  pr(8,2) , pp ( 1 6 ) 
equivalence  (p,pr) 

call  xfer(pp,pr, 16) 
anorm=0 . 0 
do  11  j=l , 8 

11  anorm=anorm+pr ( j , 1 )**2 
ortho=0,0 

do  12  j  =  1 , 4 

12  ortho=ortho+dconjg(p(j,l))*p(j,2) 
ortho=ortho/anorm 

do  13  j  =  1 , 4 

13  p(j,2)=p(j,2)-ortho*p(j,l) 
bnorm=0 . 0 

do  14  j  =  l  ,8 

14  bnorm=bnorm+pr( j ,2)**2 
anorm=l . 0/dsqrt ( anorm) 
bnorm=l .  0/dsqrt(bnorm) 
do  15  j  =  1 , 8 
pr(j,l)=pr(j,l)*anorm 

15  pr(j,2)=pr(j,2)*bnorm 
call  xfer(pr,pp, 16) 

i f ( i f 1 ag  .ne.  0)  return 

o  sum=o  sum+ortho*aprod/bprod 

aprod=aprod*anorm 

bprod=bprod*bnorm 

return 


entry  wf  stor 
1 evl =1 evl +1 

call  xfer(p  save,p  etc ( 1 , 1 evl ) , 27 ) 
o  sum=0.0 
a  prod  =  l .  0 
b  prod=l . 0 
return 

entry  wf  load 

call  xfer(p  etc ( 1 , 1 evl ) , p  save, 27) 

1 evl =1 evl - 1 

return 

end 


entry  wf  stor 


entry  wf  load 


subroutine  wf  step (p , dpdh , ht , del  h , i fl  ag) 

wf  step  increments  the  solution  of  p  from  ht  to  ht-delh, 
using  runge-kutta  integration 

i f 1  a g = 0  one  large  step,  theta 
i f 1 ag  = 1  first  small  step,  theta 
i f 1 ag  =  2  second  small  step,  theta 
iflag=3  first  small  step,  theta-dtheta 
i f 1 ag = 4  second  small  step,  theta-dtheta 

implicit  real*8  (a-h,o-z) 
common/wf  con/omega , wave  nr 
common/wf  fl ag/precsn , i so , idbg 
common/t  mtx/t(18) 
common/tm  mtx/tm(18) 
dimension  p( 16) , dpdh ( 16) , pO ( 16 ) , 

$  hdelp0(16) , del  pi ( 1 6) , del p2 ( 16) 

dimension  t  savel(18),t  save2(18), 

$  tm  savl(18),tm  sav2(18),tm  sav3(18),tm  sav4(18) 

htO=ht 

delh  k=delh*wave  nr 
hdelh  k=delh  k*0.5 
do  11  j-1,16 
P0(j)=p(j) 

hdel pO ( j )=-dpdh ( j ) *hdel h  k 
11  p(j)=pO(j)+hdelpO(j) 


ht=ht0-0 . 5*del h 
if (i flag  .le.  2) 
i f ( i f 1 ag  .eq.  0) 
i f ( i fl ag  . eq .  0) 
i f ( i f 1 ag  .eq.  1) 
i f ( i f 1 ag  . eq .  2 ) 
i f ( i fl ag  . eq .  3 ) 
i f ( i fl ag  . eq .  4) 
i f ( i fl ag  . ge  .  3) 


t  mtrx(ht) 
xfer(t,t  savel,18) 
xfer(tm,tm  sav2, 18) 
xfer (tm, tm  savl , 18) 
xfer (tm, tm  sav3 ,18) 
xfer(tm  savl , tm, 18) 
xfer (tm  sav3 , tm, 18) 
ti  mtrx 


call  p  der i v (p , dpdh ) 
do  12  j  =  1 , 16 

del  pi  ( j )  =  -dpdh ( j )*del h  k 

12  p(j)=p0(j)+0.5*delpl(j) 

call  p  deriv(p,dpdh) 
do  13  j  =  l ,  16 

delp2(j)=-dpdh(j)*delh  k 

13  p( j ) =p0 ( j ) +del p 2 ( j ) 


ht  =  htO  -del h 
i f ( i f 1 ag  .eq. 
if (i flag  .eq. 
i f ( i f 1  a g  . eq  . 
i f ( i f 1 ag  . eq . 
i f ( i f 1 ag  . eq . 
if(iflag  .eq. 
i f ( i fl ag  . eq  . 
i f ( i fl ag  . ge . 


0)  call  t  mtrx(ht) 

0)  call  xfer(t,t  save2,18) 

1)  call  xf er ( t  savel , t , 1 8) 

2)  call  xfer(t  save2,t,18) 
0)  call  xfer(tm,tm  sav4,18) 

3)  call  xfer(tm  sav2,tm,18) 

4)  call  xfer(tm  sav4,tm,18) 
3)  call  ti  mtrx 


subroutine  xfer(a,b,n) 

transfer  array  a  into  array  b. 

real*8  a,b 
dimension  a(l) ,b(l) 

do  11  j=l , n 
H  b(j)=a(j) 
return 
end 


A  A.N, 


subroutine  wf  dens  (ht,  en,  coll) 

wf  dens  computes  the  ion  density  and  collision  frequency  for  each 
specie  by  logarithmic  interpolation  of  the  corresponding  profiles, 
profile  values  are  interpolated  between  entries  mht  and  mht+1 
(Iht  and  lht+1). 

implicit  real*8  (a-h,o-z) 

common/wf  prof/enht(101) ,enlog(101 ,5) ,coll ht(25) ,col 1 f r ( 25 , 5 ) , 

$  nrlht,lht,nrmht,mht,charge(5),ratiom(5),nrspec 

dimension  en(5),  coll (5),  dele(5),  delc(5) 


if(lht  .eq.  0)  then 
1 save=0 
lht-lht+1 
if(ht  .It.  enht(lht))  go  to  10 
mht=0 
msave=0 
mht=mht+ 1 

i f (ht  .It.  collht{mht 
end  i  f 

i f ( 1 h t  .ne.  lsave)  then 


collht(mht))  go  to  20 


i f ( 1  h t  .ge.  nrlht)  1 ht=nrl ht- 1 
del h=enht(lht+l ) -enht(lht) 
do  150  k=l,nrspec 

150  dele(k)=(enlog (lht+1, k)-enlog(lht,k))/delh 
1 save  =  l  ht 
end  i  f 

if(mht  .ne.  msave)  then 

if(mht  .ge.  nrmht)  mht=nrmht-l 
delh=collht(mht+l)-collht(mht) 
do  250  k=l,nrspec 

250  delc(k)=(collfr(mht+l,k)-collfr(mht,k))/delh 
msave=mht 
end  i  f 

dh=ht-enht(l ht) 
dc  =  ht- col  1 ht (mht ) 
do  500  k=l,nrspec 

en(k)=dexp(enlog(lht,k)+dh*dele(k)) 

500  col  1 (k)=dexp (coll fr (mht , k) +dc*del c (k) ) 
return 
end 


A-19 


subroutine  wf  init(p) 


wf  init  computes  the  initial  p  matrix,  i.e.,  the  initial  conditions 
for  the  integration  dp/dz=- i k*t*p . 

implicit  rea1*8  (a-h,o-z) 

common/ t  mtx/tll,t31,t41,tl2,t32,t42,tl4,t34,t44 
common/wf  fl ag/precsn , i so , idbg 

complex*16  b3 , b2 , bl , bO , i , q  temp , det , sqroot , q (4 ) , zp , 

$  tll,t31,t41,tl2,t32,t42,tl4,t34,t44 

dimension  p ( 8 , 2 ) 
data  i/ ( 0 . dO , 1 . dO)/ 

i  f ( i s  o  .ne.  0)  go  to  50 

b3=- (til +t 44 ) 

b2  =  t 1 1 *t44 - 1 1 4*t41 - 1 3  2 

bl=-(-t32*(tll+t44)+tl2*t31+t34*t42) 

b0=-tll*(t32*t44-t34*t42) 

S  +tl2*(t31*t44-t34*t41) 

$  -tl4*(t31*t42-t32*t41) 

call  quartc ( b3 , b2 , bl , bO , q) 

qi  min=-i*q(l) 
j  1  =  1 

do  22  j=2 , 4 
q  imag=-i*q(j) 

i  f (q  imag  .gt.  qi  min)  go  to  22 
qi  min=q  imag 

j  1  =  j 

22  continue 

qr  max=q ( 1 ) 
j  2  =  1 

do  23  j  =  2 , 4 
q  real =q ( j ) 

if(q  real  .It.  qr  max)  go  to  23 
qr  max=q  real 
j2  =  j 

23  continue 

i f ( j 1  . eq .  j2)  go  to  80 

q  temp  =  q ( j 2 ) 
q( 1 )-q( jl) 
q(2)=q  temp 

do  31  j= 1 , 2 

det=(tll-q(j))*(t44-q(j))-tl4*t41 
Zp=(tl2*q(j)-(tl2*t44-tl4*t42))/det 
p(l, j)=dreal (zp) 

P ( 2 , j ) =d i mag ( zp ) 
p ( 3 , j )  =  1  .0 
P (4 , j ) =0 . 0 

zp  =  q ( J  ) 

p(5, j)=dreal (zp) 
p(6 , j) =d imag (zp) 

zp=(t42*q(j)+(tl2*t41-tll*t42))/det 
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P ( 7 , j ) =dreal (zp) 

31  p (8, j ) =d i mag ( zp) 

40  if(idbg  .It.  2)  return 
print  902, q 
print  901 , p 
return 

50  b 1 = ( t 1 l+t44)*0 . 5 
b0=tll*t44-tl4*t41 
sqroot=cdsqrt(bl**2-b0) 
q( 1 ) =bl+sqroot 
q ( 4 ) =b 1  - sqroot 
sqroot=cdsqrt(t32) 
q(2)=+sqroot 
q(3) =- sqroot 

ql  test=q( 1 ) +i *q { 1 ) 

q4  test=q(4)+i*q(4) 

i f ( q 4  test  .gt.  ql  test)  q(l)=q(4) 

q2  test=q(2)+i*q(2) 

q3  test=q(3)+i*q(3) 

if(q3  test  .gt.  q2  test)  q(2)=q(3) 

zp=tl4 

p ( 1 » 1 ) =dreal (zp) 

P  ( 2 , 1 ) =d imag ( zp ) 

P (3 , 1 ) =0 . 0 

P  ( 4 , 1 ) =0 . 0 

P ( 5 , 1 ) =0 . 0 

p(6 , 1 ) =0 . 0 

zp=- (tll-q(l)) 

P ( 7 , l)=dreal (zp) 

P (8, j ) =d imag ( zp ) 

P  ( 1 >  2 ) =0 . 0 
P ( 2 , 2 ) =0 . 0 

P ( 3 , 2 ) =0 . 0 
p ( 4 , 2 ) =0 . 0 
zp=q(2) 

p(5, 2)=dreal (zp) 
p(6,2)=dimag(zp) 
p  ( 7 , 2 ) =0 . 0 

p (8 , 2 ) =0 . 0 

go  to  40 

80  print  900, q 
stop 

900  format (’ Oerror  in  wf  init,  q  values  do  not  sort’/ 

$  ’  q=’ ,4(lpel5.5,lpel3.5)) 

901  format(’0p  values=’/ 

$  4(lpel5.5,lpel3.5)/4(lpel5.5, 1  pel  3. 5) ) 

902  format (’ Oi n i t i al  values  from  wf  init  at  ht=topht’/ 

S  ’  q= ’ , 2 ( 1 pe 1 5 . 5 , 1 pe 13 . 5 ) ) 

end 
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subroutine  quartc  (fourb3,  s i x b 2 ,  fourbl,  bO,  q) 
quartc  finds  the  roots  of  a  quartic  polynomial  from  the  closed  form, 
implicit  real*8  (a-h,o-z) 

comp! ex* 16  b3,b2,bl,b0,q,fourb3,sixb2,fourbl,b3sq,h,i,g,hprime, 

$  gprime, sqroot ,ppl us , p , 1 ogp, cbertO,cbert 1 , cbert2 , omega  1 , 

$  omega2,rootp,rootq,rootr 

real*8  mgpl us , mgmnus 
dimension  q ( 4 )  ,pri (2) 
equivalence  (p,pri) 

data  omegal/ ( - . 5d0 ,  . 86602 54038d0 ) / , omega2/ ( - . 5d0 , - . 8660254038d0)/ 
data  precsn/1 . d - 1 0/ 

b3=fourb3*0 . 25 
b2=s i xb2/6 . 0 
bl=fourbl*0 . 25 
b3sq=b3**2 
h=b2-b3sq 

i=b0-4.0*b3*bl+3.0*b2**2 
g=bl+b3*(-3.0*b2+2.0*b3sq) 
hprime=- i/12 .0 

gprime=-g**2/4.0-h*(h**2+3.0*hprime) 
sqroot=cdsqrt(gprime**2+4.0*hprime**3) 
p=(-gprime+sqroot)*0.5 
mgplus=dabs(pri(l))+dabs(pri(2)) 
ppl us=p 

p=(-gprime-sqroot)*0.5 
mgmnus=dabs(pri ( 1 ) )+dabs (pri (2) ) 
if(mgplus  .gt.  mgmnus)  p=pplus 
1 ogp=cdl og (p ) 
cbert0=cdexp(logp/3.0) 
cbertl=omegal*cbertO 
cbert2=omega2*cbert0 
rootp=cdsqrt (cbertO- hpr i me/cbertO- h ) 
rootq  =  cdsqrt(cbertl - hpr i me/cbert 1  - h ) 
rootr=cdsqrt(cbert2-hpri me/cbert 2  -  h ) 
if(cdabs(g)  .le.  1.0d-20)  go  to  5 
sign=-rootp*rootq*rootr*2.0/g 
if(sign  .It.  0.0)  rootr=-rootr 

q(l)=+rootp+rootq+rootr-b3 

q(2)=+rootp-rootq-rootr-b3 

q(3)=-rootp+rootq-rootr-b3 

q(4)=-rootp-rootq+rootr-b3 

do  20  n= 1 , 4 
i  t  e  r  =  0 

rootp=q(n)**4+fourb3*q(n)**3+sixb2*q(n)**2+fourbl*q(n)+b0 
rootq=4.0*q(n)**3+3.0*fourb3*q(n)**2+2.0*sixb2*q(n)+fourbl 
rootr=rootp/rootq 
q(n)=q(n) -rootr 

i f (cdabs (rootr )  .It.  precsn)  go  to  20 
i  ter=i ter+1 

i f ( i ter  .It.  10)  go  to  10 
print  900,  iter,  q(n) 
continue 


subroutine  rbars(c,s,rbarll,rbar22,ey,hy) 
implicit  real*8  (a-h,o-z) 

common/wf  inpt/theta,freq,azim,codip,magfld,coefnu(5),expnu(5) 
topht,lwstht,delht,h, alpha, sigma, epsr,nelect 
common/wf  con/omega, k 
common/ey  grnd/eyg,hyg 

comp 1  ex* 16  theta , i , ngsq,c,s,ssq,sq root, rtiort, ike, 
pO , hlO , h  2  0 , hlprmO , h2prm0 ,caphl0,caph20, 
pd,hld,h2d, hi prmd , h2prmd ,caphld,caph2d, 
pz,hlz,h2z,hlprmz,h2prmz, 
a  1st , a2nd , a3rd, a4th , al , a2 , a3 , a4 , 
exd,exdsq,exz,exzsq, 
rbarl 1 , rbar22 , z 1 , z  2 , 
denl2,den34, 
eyg,hyg, 

ex , ey , ez , hx , hy , hz 

real *8  k,kvraot,kvratt,nOsq,ndsq,nzsq,magfld,lwstht 
equivalence  (pz,pd),(hlz,hld),(h2z,h2d), ( h 1 prmz , h 1 prmd ) , 

( h2prmz , h2prmd ) ,(exd,exz),(exdsq,exzsq) 

data  i / ( 0 . OdO ,  1 . OdO )/ 
data  tstthm/lO.dO/ 
data  epsl n0/8.85434d- 12/ 

eyg=l  .0 
hyg=  1 . 0 
al t=l wstht 
d=l wstht 
ssq=s*s 

ngsq=dcmpl x(epsr , -  sigma/ (omega*epsl nO) ) 

sqroot=cdsqrt(ngsq-ssq) 

thtim=i*theta 

if(thtim  .gt.  tstthm)  go  to  10 

kvraot=dexp(dlog(k/alpha)/3.0) 

kvratt=kvraot**2 

avrkot=l . 0/kvraot 

avrktt=avrkot**2*0.5 

nOsq=l .0-alpha*h 

rtiort=nOsq/ngsq*sqroot 

pO=kvratt*(nOsq-ssq) 

call  mdhnkl  (pO , hi 0 , h 2 0 , h 1 prmO , h2prm0 , theta rb  1’) 

caphlO=hlprmO+avrktt*hlO 

caph20=h2prm0+avrktt*h20 

alst=caph20-i*rtiort*kvraot*h20 

a2nd=caphl0-i*rtiort*kvraot*hl0 

a3rd=h2prm0-i*kvraot*sqroot*h20 

a4th=hlprm0-i*kvraot*sqroot*hl0 

denl2=h20*a2nd-hl0*alst 

den34=h20*a4th-hl0*a3rd 

if(d  .eq.  0.0)  go  to  10 

ndsq=l .0-alpha*(h-d) 
pd=kvratt*(ndsq-ssq) 

call  mdhnkl  (pd , hid , h 2 d , hi prmd , h2prmd , theta rb  2’) 
caph 1 d  =  h lprmd+avrktt*h Id 
caph2d=h2prmd+avrktt*h2d 


al=c*ndsq*(h2d*a2nd-hld*alst) 

a2=i*avrkot*(caphld*alst-caph?d*a2nd) 

a3=i *avrkot* (h2prmd*a4th -h lprmd*a3rd ) 

a4=c*{h2d*a4th-hld*a3rd) 

rbarll=(al-a2)/(al+a2) 

rbar22=(a3+a4)/(a4-a3) 

expon=dexp(0.5*alpha*alt) 

hy=(h2z*a2nd-hlz*alst)*expon/den 1 2  *  hyg 

ey=(h2z*a4th-hlz*a3rd)/den34*eyg 

retu  rn 

flat  earth 
10  ikc=i*k*c 

exd=cdexp ( - i kc*d ) 
exdsq=exd**2 

zl=(ngsq*c-sqroot)/(ngsq*c+sqroot) 

z2=(c-sqroot)/(c+sqroot) 

rbar 1 1 =z l*exdsq 

rbar22=z2*exdsq 

hy=( 1 . 0+zl*exzsq)/ ( 1 . 0+z 1 )/exz*hyg 
ey=(1.0+z2*exzsq)/(l . 0+z 2 )/exz*eyg 
return 

entry  wf  htgn 

entry  wf  htgn ( ht , ex , ey , ez , hx , hy , hz , c , s ) 
al t=ht 

i f ( t h t i m  .gt.  tstthm)  go  to  50 
nzsq  =  l . 0  -  a  1 pha*( h- al t ) 
pz=kvratt*(nzsq-ssq) 

call  mdhnkl  ( pz , h lz , n2z , h lprmz , h2prmz , t heta  , ’htgn’ ) 

expon=dexp(0.5*alpha*alt) 

hy  =  ( h2z*a2nd -h 1 z  *  a  1 st ) *expon/den 1 2  *  hy g 

ey^(h2z*a4th-hlz*a3rd)/den34*eyg 

ex=i*avrkot*( ( h2prmz*a2nd - h 1 prmz*a 1 s t ) / 

S  denl2*hyg*expon+avrktt*hy)/nzsq 
ez=- s/nzsq*hy 
hz=s*ey 

hx=avrkot/i*(h2prmz*a4th-hlprmz*a3rd)/den34*eyg 

return 

flat  earth 

50  exz  =  cdexp ( - i kc*al t ) 
exzsq=exz**2 

hy=(1.0+zl*exzsq)/(l . 0+zi)/exz*hyg 
ey= ( 1 . 0+z2*exzsq)/ ( 1 . 0+z2)/exz*eyg 
ex=-c*(1.0-zl*exzsq)/(1.0+zl)/exz*hyg 
ez=- s*hy 
hz=s*ey 

hx=c*(1.0-z2*exzsq)/(1.0+z2)/exz*eyg 

return 

end 
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subroutine  r  mtrx(p,cosn,r) 

r  mtrx  computes  reflection  coefficient  matrix  from  p  matrix 
and  returns  it  in  r. 

implicit  real*8  (a-h,o-z) 
complex*16  zp ( 4 , 2) , cosn , r ( 2 , 2 } , 

S  g  1 2 ,  g 1 3 ,  gl4,  g23,  g24,  g34, 

S  dOO,  dll,  d 2 2 ,  d 1 2 ,  d21 

dimension  p ( 8 , 2 ) 

do  2  j  =  1 , 2 

zp ( I » j ) =dcmpl x (p ( 1 , j ) , p(2 , j ) ) 
zp ( 2 , j )=dcmplx(p(3, j) ,p(4,  j) ) 
zp(3,j )=dcmplx(p(5,j) ,  p  ( 6 ,  j  ) ) 
zp(4,j)=dcmplx(p(7,j),p(8,j)) 

gl2=zp(l,l)*zp(2,2)-zp(l,2)*zp(2,l) 
gl3=zp(l,l)*zp(3,2)-zp(l,2)*zp(3,l) 
g  1  4  =  z p ( 1  ,  l)*zp(4,2)-zp(l , 2 ) * z p ( 4 , 1 ) 
g23=zp(2,l)*zp(3,2)-zp(2,2)*zp(3,l) 
g24=zp(2,lj*zp(4,2)-zp(2,2)*zp(4,l) 
g34=zp(3,l)*zp(4,2)-zp(3,2)*zp(4,l) 

d00=-gl3+cosn*(  g34-gl2+cosn*g24) 
d  1 1  =  gl3  +  cosn*(  a34+gi2+cosn*g24) 
d  2  2  =  gl3  +  cosn*(-g34-gl2  +  cosn*g24) 
d  1 2  =  2 . 0*cosn*g 14 
d  2 1  =  2 . 0*cosn*g23 

r ( 1 , 1 ) =d 1 1/dOO 
r ( 2 , 2) =d22/d00 
r ( 1 , 2 ) =d 1 2/d00 
r ( 2 , 1 ) =d2 1/dOO 
return 
end 
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subroutine  t  mtrx(ht) 

t  mtrx  computes  m-  the  susceptibility  tensor  and 

t-  the  coefficient  matrix  of  dp/dz=- i k*t*p . 
note  that  on  call  to  entry  init  t.  various  ionospheric 
constants  are  computed. 

implicit  real*8  (a-h,o-z) 

common/wf  fl ag/precsn , i so, idbg 

common/m  mtx/ml 1 ,m21 ,m31 ,ml2 ,m22 ,  m32 ,ml3 ,  m23 ,m33 

common/ t  mtx/tll ,t31 ,t41 ,tl2,t32,t42,tl4,t34,t44 

common/tm  mtx/tml 1 , tm3 1 , tm41 , tml 2 , tm32 , tm42 , tml 4 , t m3 4 , tm44 

common/c  s/c , s , c i , s i 

common/wf  inpt/theta,freq,azmuth,codip,magfld,ceffnu(5),expnu(5), 

S  topht,lw$tht,delht,h, alpha, sigma, epslon,nglect 

common/ i trat/dtheta,dlub(2),maxitr, itr 
common/wf  con/oniega  ,  wave  nr 

common/wf  prof/enht(101),enlog(101,5),collht(25),collfr(25,5), 

$  nr1ht,lht,nrmht,mht,charge(5),ratiom(5),nrspec 

real *8  magfld,lwstht,lsqysq,msqysq,nsqysq,lmysq,lnysq,mnysq, 

S  nu  ,  1 y , my , ny 

complex*16  m(3 , 3) , 

$  mll,m21,m31,ml2,m22,m32,mi3,m23,m33, 

S  tll,t31,t41,tl2,t32,t42,tl4,t34,t44, 

S  tml 1 , tm31 , tm41 , tml 2 , tm32 , tm42 , tml 4 , tm34 , tm44 , 

S  c , s , c i , s i , esq , ssq  ,  esq i , ssqi , 

S  theta , dtheta , 

S  d,ml3d,m23d, 

$  u,usq,dd, i , iud,ta,tb 

dimension  y ( 5) ,ysq ( 5 ) , ly ( 5) , my ( 5) , ny (5) , coef  en ( 5 ) , en ( 5 ) , nu ( 5 ) , 

$  1 mysq ( 5 ) , 1 ny sq ( 5 ) ,mnysq(5) ,lsqysq(5),msqysq(5),nsqysq(5) 

equivalence  (mll,m) 

data  pi/3. 1 4 1 592653d0/ 

data  two pi/6. 2831853071795 9d0/ 

data  dtr/0.01745329252d0/ 

data  coeffx/3. 1823 57d03/ , coef fy/ 1 . 758796d 11/ 
data  i / ( 0 . dO ,  1 . dO ) / 
data  vel  1  1/2 .997 928d0 5/ 

calculate  the  matrix  m. 
m ( 1 ,  1 ) =0 . 0 
m( 1 , 2) =0 . 0 
m( 1 , 3) =0 . 0 
m( 2 , 1 ) =0 . 0 
m( 2 , 2 ) =0 . 0 
m( 2 , 3 ) =0 . 0 
m( 3 , 1 ) =0 . 0 
m(3 , 2) =0 . 0 
m(  3 , 3) =0 . 0 

call  wf  dens  (ht,  en,  nu) 
nf 1 ag=0 

do  20  k= 1 , nrspec 

add  in  the  contributions  to  the  susceptibility  tensor  m  for  each 
specie  in  the  ionosphere. 

if(en(k)  .It.  l.0d-3)  go  to  20 
n  f  1  a  g  =  1 
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x=coef  en(k)*en(k) 
if(nglect  .ne.  0)  x=-x 
z=nu(k)*ov  omga 
if(nglect  .ne.  0)  z=-z 
u  =  l .  0-  i*z 
usq=u*u 

dd=-x/(u*(usq-ysa(k) )  ) 
iud=(z+i )*dd 
ta=usq*dd 
m( 1 , 1 ) =m( 1 , 1 )+ta 
m( 2 , 2 ) =m ( 2 , 2 )+ta 
m(3,3)=m(3,3)+ta 
m(2,2)=m(2,2)-msqysq(k)*dd 
ta=my(k)*iud 
tb=l nysq(k)*dd 
m( 1 , 3 ) =m( 1 , 3 )  +  ta-  tb 
m( 3 , 1 ) -m(3 , 1 ) - ta- tb 
i f ( i so  . ne .  0)  go  to  20 
m( 1 , 1 ) =m( 1 , 1 ) -1 sqysq (k) *dd 
m(3,3)=m(3,3)-nsqysq(k)*dd 
ta=ny(k)*iud 
t  b  =  1  my sq ( k) *dd 
m ( 2 , 1 ) =m ( 2 , 1 ) +ta - tb 
m( 1 , 2 ) =m ( 1 , 2 ) - ta - tb 
ta=ly(k)*iud 
tb=mnysq(k)*dd 
m( 3 , 2 ) =m (3 , 2 )  +  ta - tb 
m(2,3)=m(2,3)-ta-tb 
20  continue 


crvtrm=al pha*(h-ht ) 
m(  1 , 1 ) =m( 1 , 1 ) -crvtrm 
m(2,2)=m(2,2)-crvtrm 
m(3,3)=m(3,3)-crvtrm 


calculate  the  matrix  t. 
d=l  .0/(1 ,0+m33) 
tm4 1 = 1 .0+mll 
tm32=m22 
tml4=d 

if(nflag  .eq.  0)  go  to  40 

ml 3d=ml3*d 

m23d=m23*d 

tm41=tm41-m31*ml3d 

tml l=m31*d 

tm44=ml3d 

i f ( i so  . ne .  0)  go  to  40 
tm32=tm32 -m32*m23d 
tm31=m31*m23d-m21 
tml2=£m32*d 
tm42=m32*ml3d-ml 2 
tm34=m23d 


40  1 4 1  =  tm4 1 
t32=csq+tm32 
1 1 4  =  1 .0-ssq*tml4 
if(nflag  .eq.  0)  go  to  70 
tl l=-s*tml 1 
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t44=-s*tm44 

i f ( i s o  .ne.  0)  return 

t31=tm3 1 

tl2=s*tml2 

t42=tm42 

t34=s*tm34 

return 

entry  ti  mtrx 
t4 1 =tm4 1 
t32=csqi +tm32 
1 1 4  =  1 .0-ssqi*tml4 

1 1 1  —  - s i * tml 1 
t44=-si*tm44 

i f ( i s o  .ne.  0)  return 
t31=tm31 

1 1 2  =  s i *tml 2 
t42=tm42 

1 3  4  =  s i *  t m3  4 
return 


entry  ti  mtrx 


entry  i n i t  t 


entry  init  t 
1  ht  =  0 
i  so  =  0 

if(magfld  .eq.  O.dO)  go  to  250 
if (dabs(codip-90.d0)  .ge.  0 . 1 5d0 )  go  to  300 
i f (dabs ( azmuth -  90. dO)  .It.  0 . 1 5d 0 )  go  to  250 
if (dabs(azmuth-270.d0)  .ge.  0.15d0)  go  to  300 

250  i s o=  1 

300  omega  =  twopi *freq*1000  .  dO 
ov  omga= 1 . 0/omega 
wavenr=omega/vel 1 t 
sindip=dsin(codip*dtr) 
drcosi=sindip*dcos(azmuth*dtr) 
drcosm=sindip*dsin(azmuth*dtr) 
drcosn=-dcos(codip*dtr) 
do  60  k=l,nrspec 

coef  en(k)=coeffx*1.0d6*charge(k)**2/(omega**2*ratiom(k)) 

y( k) =coeffy*charge(k)*magfld/(omega*rati om(k) ) 

ysq ( k) =y ( k) **2 

1y(k)=drcosl*y(k) 

my(k)=drcosm*y(k) 

ny(k)=drcosn*y (k) 

lsqysq(k)=drcosl**2*ysq(k) 

msqysq ( k) =drcosm**2*y sq { k) 

nsqysq(k)=drcosn*llr2*ysq{k) 

lmysq( k) =drcos1 *drcosm*ysq ( k ) 

lnysq(k)=drcosl*drcosn*ysq(k) 

mnysq( k) =drcosm*drcosn*y sq ( k ) 

60  continue 

c=cdcos ( theta*dtr) 
s=cds i n ( theta*dtr ) 
csq=c**2 
ssq=s**2 

ci-cdcos((theta-dtheta)*dtr) 
si=cdsin((t'neta-dtheta)*dtr) 
csqi =c i **2 
ssqi=si**2 


I 


% 

I 

1 


subroutine  wf  bndy(b) 

wf  bndy  computes  the  vector  b,  which  determines  how  to  combine 
the  solution  vectors  in  order  to  satisfy  the  boundary  conditions, 
this  routine  is  valid  only  for  eigenangles  of  the  modal  equation 
and  is  used  to  compute  height  gain  functions. 

implicit  real*8  (a-h,o-z) 
common/wf  flag/precsn, iso, idbg 
common/p  mtx/p (4 , 2 ) , pi ( 1 6 ) 
common/ey  grnd/eyg,hyg 

complex* 16  p , b (2 ) , r (2 , 2 ) , f , rbar 1 1 , r bar 2 2 , nurmf , denmf , numa , dena , a , 
S  exl ,ex2,eyl ,ey2 ,hxl ,hx2, hyl , hy2 ,ey , hy , eyg , hyg , c , s , fof r 

S  hysum , abparl , abperp 

dimension  pp ( 1 6 ) 
data  dt r/0 . 0 1 7 4 5 3 2 9 2 5 2 dO/ 

i f ( i dbg  .gt.  1)  print  902,  p 

ex  1  =  p ( 1 , 1) 

ex2  =  p ( 1 , 2 ) 

ey  1  =  -  p  { 2 , 1 ) 

ey  2  =  -  p  ( 2  ,  2  ) 

hx 1  =  p ( 3 ,  1 ) 

hx2  =  p ( 3 , 2  ) 

hyl=p(4,l) 

hy2=p ( 4 , 2 ) 

i f ( i so  . ne .  0)  go  to  500 

compute  b,  non  -  i sotrop  i  c  case  (from  polarization  ey/hy,  see  budden) 
nurmf=r(2, 1)*(1 . 0+rbar22 )*rbar  1 1 
denmf = (1.0+rbarll)*(1.0-rbar22*r(2,2)) 
f  of  r=nurmf/denmf 
numa  =  - (ey 1  - fof r*hy 1 ) 
dena  =  ey2  -  f of r*hy2 
a=numa/dena 
hysum=hy l+a*hy2 
b( 1 ) =1 . 0/hysum*hy 
b(2)=a/hysum*hy 
eyg= (ey 1 *b ( 1 ) +ey2*b ( 2 ) )/ey 
go  to  820 

compute  b,  isotropic  case  (choose  correctly  polarized  solution). 

500  abparl =1 .0-rbarl  l*r( 1 , 1) 
abperp=1.0-rbar22*r(2,2) 
temp  a=cdabs (abperp) 
temp  b=cdabs (abparl ) 
temp=temp  a/temp  b 
if(temp  .It.  l.OdO)  go  to  600 
b(  1 )  =  1  .  0/hy 1 *hy 
b ( 2 ) =0 . 0 
eyg  =  0 . 0 
go  to  700 

600  b( 1 ) =0 . 0 

b ( 2 ) = 1 . 0/ey2*ey 
hyg  =  0 . 0 

700  if(temp  .It.  lO.OdO)  go  to  800 

if(temp  .gt.  O.ldO)  go  to  800 

if(temp  b  .gt.  O.ldO)  go  to  800 
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if(temp  a  .gt.  O.ldO)  go  to  800 
i f ( i dbg  .eq.  0)  go  to  820 
800  print  900,  abparl ,  abperp 

820  continue 

i f ( i dbg , ge . 2 )  print  905,  b 
return 


entry  f  fct 


entry  f  fct(pp,c,s,f ) 

cal  1  r  mtrx(pp ,c , r) 

call  rbars(c,s,rbarll,rbar22,ey,hy) 

if(idbg  .gt.  1)  print  904,  r,  rbarll,  rbar22 

compute  modal  eqn.  value 
a= ( 1 . 0 -r ( 1 , 1 )*rbar 1 1 ) 
a=a*(1.0-r(2,2)*rbar22) 
f=a- r ( 1 , 2 ) *r (2 , 1 )*rbar 1 l*rbar22 
print  90 1 ,  f 
return 

900  f ormat ( ’ Opol ar i zat i on  values’/ 

S  ’  abparl  =  ’ , ( 1  pel 5 . 5 , 1 pel3 . 5)  ,  ’  abperp= ’ , ( ] pel  5 . 5 , 1  pel 3 . 5 ) ) 

901  format (’ Omodal  eqn  val ue  =  ’ , ( lpel 5 . 5 , 1  pel 3  .  5) ) 

902  format(’0p  values  on  entry  to  wf  bndy’/ 

$  4(lpel5.5,lpel3.5)/4(lpel5.5,lpel3.5)) 

904  format (’ 0  r= ’ , 2 ( 1 pe 1 5 . 5 , lpel 3 . 5)/6x , 2 ( 1  pel  5 . 5 , 1 pe 1 3 . 5 ) / 

S  ’Orbar=’ ,2(lpel5.5,lpel3.5)) 

905  format (’ Osol ut i on  combination  factors’/ 

S  ’  bl=’,lpei5.5,lpel3.5,’  b2= ’ , 1 pe 1 5 . 5 , 1 pe 13 . 5 ) 

end 
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subroutine  itrate 


itrate  is  the  control  routine  for  finding  an  angle,  theta,  which 
satisfies  the  modal  equation. 

implicit  real*8  (a-h,o-z) 

common/wf  inpt/theta,freq,azmuth,codip,magfld,ceffnu(5),expnu(5) 
$  topht,lwstht,delht,h, alpha, sigma, epslon,nglect 

common/i trat/dtheta , bnd  rl,bnd  im,maxitr, itr 
common/p  mtx/p ( 1 6) , pi ( 16 ) 
common/cs/c , s , ci , si 
real*8  magfld,lwstht 

complex*16  theta , c , s , c i , si ,dtheta , i , f,  fO , dfdt , del  t 
data  i/(0.d0,l.d0)/ 

nr  iter=0 

i f ( i tr  . eq .  1 )  then 

11  call  wf  i ntg (topht , 1 wstht ,  del ht ,  1 ) 
cal  1  f  f c t ( p i ,ci ,si , f 0 ) 
cal  1  f  fct (p , c , s , f ) 
dfdt=(f-fO)/dtheta 
del  t=-f/dfdt 
theta=theta+del  t 
print  900, theta 
nr  iter=nr  iter+1 
if(nr  iter  .gt.  maxitr)  then 
print  *,’max  itr  exceeded’ 
stop 
end  if 

del  rl=del  t 

if(dabs(del  rl )  .gt.  bnd  rl )  go  to  11 
del  im=-i*del  t 

if(dabs(del  im)  .gt.  bnd  im)  go  to  11 
end  i  f 

call  wf  intg(t op h t , 1 wstht, del ht,0) 

cal  1  f  fct (p , c , s , f ) 

return 


)  format(’0new  theta  =  ’ , 2 f 1 0 . 3 ) 


subroutine  p  deri v(p,dpdh) 


c 

c  p  deriv  computes  the  height  derivatives  of  the  field  vectors, 
c  p,  according  to  clemmow  and  heading  (1954). 
c  equation  is  dp/dz=- i k*t*p . 

c  multiplication  by  -i  is  performed  by  operating  on  real  and  imag  parts, 
c  multiplication  by  k  is  performed  in  routine  wf  step, 
c 

common/t  mtx/tll,t31,t41,tl2,t32,t42,tl4,t34,t44 
complex*16  zp(4), deriv, 

S  til , t 3 1 ,t41 ,t!2,t32,t42,t!4,t34,t44 

real*8  p (8 , 2 ) , dpdh (8 , 2 ) , part ( 2  ) 
equivalence  (deriv, part) 

do  11  j-1,2 

zp(l)=dcmplx(p(l,j),p(2,j)) 

zp ( 2 ) =dcmpl x ( p ( 3 , j ) ,p(4,  j) ) 

zp(3)=dcmplx(p(5,j) ,p(6,j)  ) 

zp(4)=dcmplx(p(7, j) ,p(8, j) ) 

deriv=tll*zp(l)+tl2*zp(2)+tl4*zp(4) 

dpdh ( 1 ,  j)=  part ( 2 ) 

dpdh ( 2 , j )  =  -  part ( 1 ) 

der i v  =  zp (3  ) 

dpdh ( 3 , j )=  part ( 2 ) 

dpdh ( 4 , j )  =  -  part ( 1 ) 

deriv=t31*zp(l)+t32*zp(2)+t34*zp(4) 

dpdh ( 5 , j ) =  part ( 2 ) 

dpdh (6 , j )=-part( 1 ) 

deriv=t41*zp(l)+t42*zp(2)+t44*zp(4) 

dpdh ( 7 , j ) =  part ( 2 ) 

11  dpdh(8,j)=-part(l) 
return 
end 
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function  cdang(arg) 
implicit  real*8  (a-h,o-z) 
complex*16  arg,mi/(0.d0, -1 .dO)/ 

argr=arg 

argi=arg*mi 

cdang=datan2(argi , argr) 

return 

end 


A  BRIEF  DESCRIPTION  OF  USAGE  OF  PROGRAM  WFPLTS 

If  the  namelist  logical  variable  SAVPLT  is  equal  to  .TRUE,  (which  is  the  default  value), 
data  are  written  to  unit  10  for  optional  further  processing  and  plotting  by  the  program  WFPLTS. 
WFPLTS  has  two  sources  of  input:  the  data  file  generated  by  WAVFLD  and  namelist  input. 
The  program  prompts  the  user  for  the  name  of  the  WAVFLD  data  file  and  then  asks  for  the 
namelist  data.  There  are  four  variables  in  the  namelist:  the  size  of  the  x-axis,  the  size  of  the 
y-axis  in  inches,  a  scaling  factor  and  a  plotting  flag.  The  plotting  flag  is  NPLOT,  an  array  of  six 
integers.  It  controls  which  plots  are  generated.  A  maximum  of  six  plots  may  be  generated. 

If  NPLOT(l)  is  non-zero  Q  and  ATTEN  are  plotted. 

If  NPLOT(2)  is  non-zero  the  magnitude  of  SX.  SY  and  SZ  are  plotted. 

If  NPLOT(3)  is  non-zero  angles  between  SX  and  SZ,  SY  and  SZ,  B  DOT  S  are  plotted. 

If  NPLOT(4)  is  non-zero  X,  Y,  and  Z  are  plotted. 

If  NPLOT(5)  is  non-zero  the  magnitude  of  EX,  EY,  and  EZ  are  plotted. 

If  NPLOT(6)  is  non-zero  the  magnitude  of  HX,  HY,  and  HZ  are  plotted. 

When  the  program  prompts  the  user  for  the  namelist  input,  default  values  for  all  variables 
are  listed,  so  only  variables  the  user  wishes  to  change  need  be  entered. 


program  to  plot  output  file  from  wavfld 

nplot  controls  groups  of  plots: 
nplot(l)  ne  0  gives  q  and  atten 

2  magnitude  of  sx,  sy  and  sz 

3  angles  between  sx  and  sz,  sy  and  sz,  b  dot  s 

4  x,  y  and  z 

5  magnitude  of  ex,  ey  and  ez 

6  magnitude  of  hx,  hy  and  hz 

dimension  pi otx ( 1 0 1 , 1 6 ) , pi oty ( 101) , u p ( 101) , x 1 (2) , y 1 (2) , u 1 (2) , 
npl  o  t  ( 6 ) , i dat e ( 3 ) , i t i me ( 2 ) , i label ( 15) , i b  c  d (20) 
character*60  pltlbl 
character*40  bed 
logical  up,  til 


namel  i  st/datum  's i/ex  .  i  /** * 

data  n  p 1 o  t  '  1 , 2 , 4  *  0  .size*, 
up  '  1 0  1  *  -  f  a  1  .  .  ■».  1  0  . 

equivalence  (  i  1  a  b>‘  1  .  ;•  1  t  1 : 


print  *  ,  *  t  n  t 
accept  1 0  0  0 . f  : 
open  (  ii n  i  t  -10,''  r 

print  *  ,  ’  f  lit  *-r 
read  datum 


begin  plot*  ; -  , 
read  (  1  0  .  err!  99 ' 


ymax=a i nt ( ht  max  10.  ♦  .99 
ym i n  =  a i n t ( ht mi  n  10.  ,  *  i  0 
ysca 1 e= ( ymax  ym i n ) 


do  2  99  rr, p  1  o t  1  .6 
i  f  (npl  ot  (rpl  of  ;•  .  •  ■  . .  0 ;  . 

call  pit.  bgn 

call  symbo 1 ( 0 . 0 , 0 . 0 , 0 . 1 . ’ wa v f 1 d‘ . 0 . . c . 
call  symbol(l-0,0.0,0.1.idate,0. .1?) 
call  symbol (2. 5,0.0,0. 1 , it ime.0. ,8) 
call  p  1  o  t  (  .5.  1  .  ,  -3) 
xp  =  0 . 0 
yp  =  -  .4 

call  symbol ( xp ,yp ,. 1 , i bed  ,0.,68) 
encode(60, 2000, pltlbl )  frq,rho,azm,cdp,$ig 
yp  =  yp-  .2 

call  symbol ( xp , y p , . I , ilabel ,0. .60) 
yp  =  yp-0 .  2 

go  to  (100, 103, 106, 110, 115, 120), mplot 


amax= - 1000 . 
do  101  i  =  1 . n  rh  t 
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plotx( j , 2)=aloglO(amaxl ( 1 .e- 10,attn-pl  otx( j ,  2) ) ) 
do  101  k=l , 2 

if(amax  .It.  plotx(j,k))  amax=pl otx( j , k) 
xmax=aint( am ax/scale x+.99)*scalex 
xmin=xmax-sizex*scalex 
xscal  e  =  s>,  - 1  ex 

call  border( si zex, xmi n , xmax, xscal e , 1 ,  si zey , ymi n ,ymax , 10 . , 1 ) 

do  102  j=l,nrht 
do  102  k=l ,  2 

i f (pi  otx ( j , k)  .It.  xmin)  pi otx  ( j , k) =xmi n 
call  symbol (xp,yp, . 1 , ’  q’,0.,2) 
xp=xp+0 . 25 

call  curve (xl ,yl , ul , 2 , -xp , -yp, 1 . , 1 . , 1 ) 
xp=xp+l .  2 

call  curve (pi otx( 1 , 1 ) , pi oty , up , nrht , xmi n ,ymi n , xscal e ,y seal e , 1 ) 
call  newpen(2) 

call  symbol (xp ,yp, . 1 , ’ at ’ ,0 .  ,2 ) 
xp=xp+0 . 25 

call  curve (xl , yl , ul , 2 , -xp , -yp , 1 . , 1 . , 4) 
xp=xp+l . 2 

call  curve(plotx(l,2),ploty,up,nrht,xmin,ymin,xscale,yscale,4) 
go  to  298 

amax=- 1000 . 
do  104  j=l,nrht 
do  104  k=3 , 5 

i f (amax  .It.  plotx(j,k))  amax=pl otx( j , k) 
xmax=aint(amax/scalex+.99)*scal ex 
xmin=xmax-si zex*scal ex 
xscal e=scal ex 

call  border(si zex, xmin, xmax, xscale,  1 ,  s i zey ,ymi n ,ymax ,  10  .  ,  1 ) 

do  105  j=l ,nrht 
do  105  k  =  3 , 5 

i f ( pi otx ( j , k)  .It.  xmin)  plotx( j , k)=xmin 
call  symbol (xp.yp, .1 ,’ sx’ ,0. ,2) 
xp=xp+0 . 25 

call  curve (xl ,yl , ul , 2, -xp, -yp, 1 .  ,  1  . ,  1 ) 
xp=xp+l . 2 

call  curve(plotx(l,3),ploty,up,nrht,xmin,ymin,xscale,yscale,l) 
call  newpen(2) 

call  symbol (xp ,yp, . 1 ,’ sy’ ,0.  ,2) 
xp=xp+0 . 25 

call  curve (xl ,yl ,ul , 2, -xp, -yp, 1 . , 1 . , 2) 
xp=xp+l . 2 

call  curve(plotx(l,4),ploty,up,nrht,xmin,ymin,xscale,yscale,2) 
call  newpen(3) 

call  symbol (xp ,yp, . 1 ,’ sz ’, 0 ., 2 ) 
xp  =  xp  +  0 . 25 

call  curve ( xl , yl , ul , 2 , -xp , -yp , 1 . , 1 . , 4) 
xp=xp+l . 2 

call  curve(plotx(l,5),ploty,up,nrht,xmin,ymin,xscale,yscale,4) 
go  to  298 

xmi n=- 180 . 
xmax=  180. 
xscal e=360 ./si zex 


call  border(sizex,xmin,xmax,xscale, 1 ,  s i zey , ymi n ,ymax , 10 . , 1 ) 
call  symbol (xp ,yp , .l.’atan(sz/sx)’ ,0. ,11) 
xp=xp+l . 15 

call  curve(xl,yl,ul,2,-xp,-yp,l.,l.,l) 
xp=xp+l . 2 

call  curve(plotx(l ,6) ,ploty,up,nrht,xmin,ymin,xscale,yscale, 1 ) 
call  newpen(2) 

call  symbol(xp,yp,.l,’atan(sz/sy)’,0.,ll) 
xp=xp+l . 1 5 

call  curve(xl , y 1 , u 1 ,2,-xp,-yp,l. ,1 . ,2) 
xp=xp+l . 2 

call  curve(plotx(l ,7) ,ploty,up,nrht,xmin,ymin,xscale,yscale,2) 
call  newpen(3) 

call  symbol (xp,yp, . 1 acos(b  .  s)’,0.,ll) 
xp=xp+l . 1 5 

call  curve(xl ,yl ,ul ,2,-xp,-yp,l.,l.,4) 
xp  =  xp+ 1 .  2 

call  curve(plotx(l,8),ploty,up,nrht,xmin,yinin,xscale,yscale,4) 
go  to  298 

amax= -  1000 . 
do  111  j  =  l  ,nrht 
do  111  k-9,10 

if(amax  .It.  plotx(j,k))  amax=pl otx ( j , k) 
xmax=aint(amax/scalex+.99)*scalex 
xmin=xmax-sizex*scalex 
xscal e=scal ex 

call  border(sizex,xmin,xmax,xscale, 1 ,  s i zey , ymi n ,ymax , 1 0 . , 1 ) 
do  112  j=l vnrht 
do  112  k=9 , 1 0 

if (plotx( j ,k)  .It.  xmin)  pi otx ( j , k ) =  xmi n 
call  symbol (xp ,yp, . 1 x 0 ., 1 ) 
xp  =  xp  +  0 . 1 5 

call  curve(xl ,yl , u 1 , 2 , -xp , -yp , 1 . ,1 . ,5) 
xp  =  xp+l .  2 

call  curve(plotx(l ,  9) ,ploty,up,nrht,xmin,ymin,xscale,yscale,5) 
call  newpen(2) 

call  symbol ( xp , yp ,. 1 z 0 ., 1 ) 
xp=xp+0 . 1 5 

call  curve ( xl , yl ,ul,2,-xp,-yp,l.,l.,3) 
xp=xp+l . 2 

call  curve(plotx(l,10),ploty,up,nrht.xmin,ymin,xscale,yscale.3) 
call  newpen(3) 

call  symbol (xp ,yp, . 1 , ’y’ ,0. , 1 ) 
xp=xp+0 . 1 5 

call  curve(xl,yl,ul,2,-xp,-yp,l.,l.,l) 

xp=(alogl0(capy)-xmin)/xscale 

i  f ( x  p  .le.  0.)  go  to  298 

i f ( x p  . g e .  sizex)  go  to  298 

call  plot(xp,  0.,3) 

call  pi ot ( xp , s i zey , 2 ) 

go  to  298 

amax  =  -  1000 . 
do  116  j=l ,nrht 
do  116  k  =  1 1 , 13 

if(amax  .It.  plotx(j,k))  amax  =  pl otx( j  ,  k) 
xmax=aint(amax/scalex+.99)*scalex 
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xmin=xmax-sizex*scalex 
xscal e=scal ex 

call  border ( s i zex , xmi n , xmax , xs cal e  ,  1 ,  sizey,ymin,ymax, 10. , 1) 
do  117  j  =  1 , nrh t 
do  117  k= 1 1 , 1 3 

if (plotx( j , k)  .It.  xmi n )  pi otx ( j , k) =xmi n 
call  symbol (xp,yp, . 1 , ’ex’ ,0. ,2) 
xp=xp+0 . 25 

call  curve(xl,yl,ul,2,-xp,-yp,l.,l.,l) 
xp=xp+1.2 

call  curve(plotx(l,ll),ploty,up,nrht,xmin,ymin,xscale,ysca!e,l) 
call  newpen(2) 

call  symbol (xp.yp, . 1 , ’ey’ ,0. ,2) 
xp=xp+0 . 25 

call  curve(xl , yl ,ul ,2,-xp,-yp,l.,l.,3) 
xp=xp+ 1 . 2 

call  curve(plotx(l,12),ploty,up,nrht,xmin,ymin,xscale,yscale,3) 
call  newpen(3) 

call  symbol ( xp ,yp ,. 1 ez 0 ., 2 ) 
xp=xp+0 . 25 

call  curve(xl ,yl ,ul ,2,-xp,-yp,l. , 1 . , 4) 
xp=xp+l . 2 

call  curve(plotx(l , 13) ,ploty ,up,nrht ,xmin,ymin,xscal  e,yscale,4) 
go  to  298 

amax  =  -  1000 . 
do  121  j=l,nrht 
do  121  k=l 4 , 16 

i f (amax  .It.  plotx(j,k))  amax=pl otx( j , k) 
xmax=aint(amax/scalex+.99)*scalex 
xmin=xmax-sizex*$calex 
xscal e=scal ex 

call  border ( s i zex , xmi n , xmax , xscal e ,  1 ,  si zey , ymi n ,ymax ,  1 0  .  ,  1 ) 
do  122  j=l,nrht 
do  122  k=14 , 16 

if(plotx(j,k)  .It.  xmin)  plotx( j , k)=xmin 
call  symbol (xp,yp, . 1 , ’hx’ ,0.  ,2) 
xp=xp+0 . 25 

call  curve(xl,yl,ul,2,-xp, -yp ,  1 . ,  1 . ,  1 ) 
xp=xp+l . 2 

call  curve(plotx(l,14),ploty,up,nrht,xmin,ymin,xscale,yscale,l) 
call  newpen(2) 

call  symbol (xp ,yp, . 1 hy’ , 0. , 2 ) 
xp=xp+0 . 25 

call  curve(xl,yl,ul,2,-xp,-yp,l.,l.,3) 
xp=xp+l . 2 

call  curve (pi otx ( 1 ,15),ploty,up,nrht,xmin,ymin,xscale,yscale,3) 
call  newpen(3) 

call  symbol (xp ,yp ,. 1 hz 0 ., 2 ) 
xp=xp+0 . 25 

call  curve (xl ,yl,ul,2,-xp,-yp,l.,l.,4) 
xp=xp  +  l .  2 

call  curve (pi otx ( 1 , 16) ,pl oty , up, nrht , xmi n,ymi n , xscal e,yscal e , 4 ) 

call  pltend 
continue 
go  to  90 


subroutine  curve(x,y,up,nrpts,xmir.,ymin,xinc,yinc,line) 
c 

c  x,y,up  must  be  dimensioned  at  least  nrpts 
c  xmin,ymin  are  x,y  origin  in  user  units 
c  xinc.yinc  are  x,y  scales  in  user  units  per  inch 
c 

c  1 i ne  =  l :  sol  id 

c  2:  long  dash 

c  3:  medium  dash 

c  4:  short  dash 

c  5:  dotted 

c  6:  short  +  long  dash 

c  7:  short  +  short  +  long  dash 

c 

logical  up,upl,up2 

dimension  ipen(10) , joc(7) , x(nrpts) ,y(nrpts),up(nrpts) 
data  i pen/3 , 2 , 3 , 2 , 3 , 2 , 2 , 2 , 2 , 2/ , joc/18,  61,  56,  54,  52,  11,  36/ 
data  delr/.l/ 
c 

i f (nrpts  . 1 e .  1 )  go  to  99 
c 

if (1 ine)  1,2,3 

1  kk=mod(l ine,7)+7 
go  to  4 

2  kk=0 

go  to  4 

3  kk=mod(l ine, 7) 

4  kk=kk+l 
jo=joc ( kk ) / 1 0 
jc=joc (kk) - 10*jo 

c 

j  =  l 

ip=2 

if(kk  .eq.  6)  i p  =  3 
dr  =  0 . 
rhol=0 . 
rho2=del r 

px  1  =  (x( 1 ) -xmi n )/xi nc 
pyl=(y(l)-ymin)/yinc 
upl-up(l) 
if(upl)  go  to  10 
c 

c  go  to  first  position  with  pen  up 
call  plot(pxl ,pyl ,3) 
c 

10  do  40  i  =  2 , nrpts 

px2=(x(i )-xmin)/xinc 
Py2=(y(i )-ymin)/yinc 
up2=up(i) 
i f ( u p 2 )  go  to  22 
if(upl)  go  to  37 
if(kk  .eq.  2)  go  to  38 
del x=px2 -pxl 
dely=py2-pyl 

rho=sqrt(delx**2+dely**2) 

rhol=rhol+rho 

if(rho2  .gt.  rhol)  go  to  38 
del x=del x*del r/rho 


de1y=dely*del r/rho 
dx  6=delx*.l 
dy  6=dely* . 1 
i f (dr  .eq.  0 . )  go  to  20 
dx=del x*dr/del r 
dy=del y*dr/del  r 
pxl=pxl+dx 
py l=py 1+dy 
go  to  21 

20  if(rho2  .gt.  rhol)  go  to  38 
pxl=pxl+delx 

pyl=py 1+dely 

21  cal  1  plot(pxl , pyl , ip) 

if(kk  .eq.  6)  call  plot(pxl+dx6,pyl+dy6,2) 
j  =  j  +  l 

i  p  =  i pen ( jo+mod( j , jc) ) 
rho2=rho2+del r 
go  to  20 

22  dr=0 . 
rho 1 =0  . 
rho2=del r 
go  to  39 

c  pen  has  been  up,  prepare  to  lower  pen 

37  call  plot(px2,py2,3) 

go  to  39 

38  call  plot(px2,py2, ip) 

dr=rho2- rhol 

39  pxl=px2 
pyl=py2 
up  1  =  up2 

40  continue 

99  return 


c 


subroutine  border(x1ng,xmin,xmax,xinc,nx,ylng,ymin,ymax,yinc,ny) 


dimension  x i nc ( nx ) , y i nc ( ny ) 
logical  fy , f x 

fx= .false. 
fy= .false. 

if(nx  .eq.  1)  fx=.true. 
if(ny  .eq.  1)  fy=.true. 
x t  =  x 1 ng-  .  1 
y t  =y 1 ng-  .  1 

xscale=xlng/(xmax-xmin) 
yscal e=yl ng/ (ymax-ymi n) 
ym=abs (ymi n ) 
yl n=  - .  4 

i  f (ym  . ge .  10 . )  yl n  =  yl n - . 1 
if(ym  .ge.  100.)  yln=yln-.l 
if (ym  .ge.  1000.)  yln=yln-.l 
if(ymin  .It.  0.)  yln=yln-.l 
ym=abs (ymax) 
yl  m= - .  4 

if(ym  .ge.  10.)  ylm=ylm-.l 
if(ym  .ge.  100.)  ylm=ylm-.l 
if(ym  .ge.  1000.)  ylm=ylm~.l 
if(ymax  .It.  0.)  ylm=ylm-.l 
xm=abs ( xmax) 
xl m=  - .  3 

if(xm  .ge.  10.)  xlm=xlm-.l 
if(xm  .ge.  100.)  xlm=xlm-.l 
i  f  (xm  .ge.  1000.)  xlm=xlm-.l 
if(xmax  .It.  0.)  xlm=xlm-.l 
if(fx)  dx=xinc(l) 
if(fy)  dy=yi nc ( 1 ) 

iy  =  l 

yl  =o . 

call  number (yl n , 0 .,. 1 ,ymin , 0  .,  1 ) 
call  plot(0. ,0 . ,3) 
if(fy)  go  to  110 

10  yp= (y i nc ( iy) -ymi n ) *y seal e 
go  to  111 

110  yl=yl+dy 
yp=yl *yscal e 

111  i f (yp  .It.  0.)  go  to  99 
if(yp  .ge.  ylng)  go  to  11 
call  plot(0. ,yp,2) 

call  plot( . 1 ,yp,2) 
call  pi ot (0 . ,yp,2) 
if(fy)  go  to  110 
iy=iy+ 1 

i f ( iy  . 1 e .  ny)  go  to  10 

11  call  p 1 o t ( 0 . , ylng, 2) 

call  number (yl m,yl ng-.l,.l,ymax,0.,l) 
call  pi ot (0  .  , y 1 ng , 3 ) 
i  x  =  1 
xl  =0 . 

if(fx)  go  to  112 

12  xp= ( xi nc ( i x) -xmi n) *xscal e 
go  to  120 

B-10 


xl =x  I  +dx 
xp  =  xl *xsca1 e 

i f (xp  . 1 t .  0 .)  go  to  99 

i  f  ( x  p  .  g  e  .  xlng)  go  to  13 
call  pi ot ( xp  ,y 1 ng , 2 ) 
call  plot(xp,yt,2) 
call  pi ot ( xp ,y 1 ng , 2) 
if(fx)  go  to  112 
ix=ix+l 

if(ix  . 1 e .  nx)  go  to  12 
call  pi  ot ( xl ng ,yl ng , 2 ) 
if(fy)  go  to  130 
i  y  =  i  y  - 1 

i f ( iy  . 1 e .  0 )  go  to  15 

yp=(yinc(iy)-ymin)*yscale 

go  to  14 

yl  =y 1 -dy 

y p  —  y 1 *yscal  e 

i  f (yp  . 1 e .  0 . )  go  to  15 

call  pi ot (xl ng ,yp, 2) 

call  pi ot ( xt ,yp , 2 ) 

call  pi  ot ( xl ng , y p , 2 ) 

if(fy)  go  to  130 

go  to  113 

call  pi ot (xl ng , 0 . , 2 ) 

call  number(xlng+xlm,-.2,.l,xmax,0.,l) 
call  pi  ot ( xl ng , 0 . , 3 ) 
if(fx)  go  to  150 
i  x  =  i  x  - 1 

i f ( i x  . 1 e .  0 )  go  to  17 

xp=(xinc(ix)-xmin)*xscale 

go  to  16 

xl =  x 1 -dx 

xp=xl *xscal e 

i f (xp  . 1 e .  0  . )  go  to  17 

cal  1  plot(xp,0.  ,2) 

call  p 1 o t ( xp  ,  .  1 , 2) 

call  pi ot ( xp , 0 . , 2 ) 

if(fx)  go  to  150 

go  to  115 

call  plot(0.,0.,2) 

call  number ( 0 2 1 , xmi n , 0 ., 1 ) 

return 

print  100, xl ng , xmi n, xmax, x i nc( 1 ) , nx,yl ng ,ymi n,ymax ,y i nc ( 1 ) , ny 
format(’0***  Error  in  BORDER:  xlng,  xmin,  xmax,  xinc(l).  nx  ='. 
S  1 p4e 1 5 . 5 , i 5/ 2 4x , ’ yl ng ,  ymin,  ymax,  yinc(l),  ny  = ’  .  1 p 4 e  1 5 . 5 

S  i 5/ ’ 0***’ ) 

call  pi  tend 
stop 
end 
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